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CHAPTER  1.   INTRODUCTION 

1.1   Historical  Introduction 

The  heuristic  approach  to  solve  mathematical  problems  in  physics 

has  attracted  attention  of  physicists  since  the  advent  of  digital  computers 

Undoubtedly,  J.  von  Neumann  was  the  first  to  propose  such  approach.   In 

1944  he  solved  a  simple  hydrodynamical  problem  numerically  by  using  the 

then  available  punch-card  equipment  at  the  Ballistic  Research  Laboratory 

which  was  developed  by  IBM  for  calculating  firing  and  bombing  tables. 

(Note  that  the  first  electronic  computer  ENIAC  was  yet  to  be  completed  at 

that  time!)   He  solved  the  one  dimensional  hydrodynamical  equation  for 

shock  waves  in  Lagrangian  coordinates: 

x    =  F' (x  )x 
tt       ay  aa 

where  x(a,t)  is  the  position  of  a  moving  fluid  particle  with  Lagrangian 

coordinate  a  at  time  t,  and  F  (x  )  gives  the  pressure  vs  specific  volume 

relation  (homentropic  equation  of  state)  .   His  idea  was  to  discretize  the 

spatial  coordinate  a  in  order  to  reduce  the  above  partial  differential 

equation  to  a  set  of  ordinary  differential  equations  of  the  form: 

m  d  x  /dt2  =  f(x    -x  )  -  f(x  -x   ,).  (1.1) 

n  n+1   n       n   n~l 

This  is  nothing  but  the  equation  of  motion  for  an  array  of  identical- 
mass  particles  connected  by  springs.   (Fig.  1.1)   He  observed  a  shock 

x  2  * 

wave  formed  in  such  mechanical  system  for  f(x)  =  (1  -  —)    ,0<x<2. 

This  was  the  first  numerical  experiment  ever  made  on  a  computing  device. 
He  further  discussed  that  the  extensive  but  well  planned  computational 
efforts  would  break  the  deadlock  encountered  by  purely  analytical  approach 


i< 

Since  there  is  no  dissipative  term  in  (1.1),  the  shock  wave  which 
von  Neumann  observed  is  essentially  a  solitary  wave  which  is  to  be 
described  in  the  next  section.   One  of  Von  Nermann's  results  is 
rep   !>;ced  i  n  Y\  g.  I .  ?  . 


|-/wM>^R^<>^rM>^H^ 


-^0?RPO^rjp-0 


Fig.  1.1.   Mass-spring  system. 


"time 


12.    3    *    5"    6    7  8    f  tO  11  12  13  11  IS  16J7  i*   If  20  U  22  23  2i  25  2Q  27  28  XJ     % 


Fig.     1.2.       Shock   wave    in    the    systerr    (1.1) 


(2) 
in  such  problems  as  turbulence  and  shock  waves.      The  next  great  step 

was  made  by  Fermi,  Pasta  and  Ulam,  who  studied  the  behavior  of  the  vibra- 
tion in  a  nonlinear  string  on  MANIAC  I  computer  at  Los  Alamos  scientific 

(3) 
laboratory  in  1953.      Their  model  was  also  a  discrete  one;  the  equation 

of  motion  is  identical  to  (1.1)  with  the  function  f  given  as: 

2 

(1)  f(x)  =  x+ax  (quadratic  nonlinearity) 

3 

(2)  f(x)    =   x+f3x  (cubic    nonlinearity) 

(3)  f(x)    =    6   x  jxj    <  d 

6  x+c   [xj  >  d     (broken  linear  function) 

where  a,  j3,  6,  ,  6,, ,  c  and  d  are  given  parameters.   Their  intention  was  to 
observe  the  mixing  of  energy  among  the  normal  modes  of  the  system  due  to  a 
weak  nonlinearity.   The  result  was,  however,  quite  contrary  to  what  they 
expected:   starting  with  the  initial  energy  given  to  the  lowest  mode,  only 
a  few  modes  were  excited,  after  which  almost  all  the  energy  was  returned 
to  the  original  mode.   This  result  cast  doubt  on  the  long-believed  hypothesis 
that  only  a  slight  nonlinearity  will  lead  a  dynamical  system  to  a  statistic- 
ally equilibrium  state.   This  numerical  experiment  by  Fermi,  Pasta  and  Ulam 
will  be  referred  to  as  FPU  throughout  this  paper. 

The  FPU  experiment  has  stimulated  extensive  research  during  the 
last  two  decades,  both  computationally  and  analytically,  among  physicists 
and  mathematicians.   The  research  directions  may  be  categorized  into 
several  areas: 

(1)  ergodic  (or  non-ergodic)  behavior  of  nonlinear  vibrating  systems, 

(2)  theory  of  the  normal  mode  coupling, 

(3)  thermal  conductivity  in  the  anharmonic  crystals 


(4)   energy  transport  properties  of  nonlinear  vibrating  systems  (especially 
in  relation  to  the  continuum  approximation) . 

Oi  course  these  areas  are  not  independent  of  each  other,  but  closely 
related.   Here,  particular  attention  should  be  paid  to  the  last  area  because 
a  special  type  of  solution  has  been  discovered  in  such  dynamical  systems; 

that  is,  the  solitary  wave  (or  the  "soliton"  as  was  named  by  Zabusky  and 

(4) 
Kruskal    ).   As  a  matter  of  fact  the  solitary  wave  has  been  known  in 

hydrodynamics  since  1844,  when  Scott  Russel  reported  a  single  elevation  of 

water  in  a  long  canal  which  propagated  without  changing  its  shape. 

Theoretical  investigations  that  followed  by  Lord  Rayleigh,     Boussinesq, 

Kortweg  and  deVries     clarified  the  nature  of  such  wave;  especially  Kortweg 

and  deVries  derived  the  equation  which  is  now  known  as   KdV  equation: 

u   +  uu   +  Bu     =  0  (1.2) 

t        X        XXX 

where  u  is  the  velocity  of  the  water.   The  maintenance  of  the  wave  is 
possible  due  to  the  balancing  of  two  different  effects:   the  nonlinearity 
that  tends  to  sharpen  the  wave  front,  and  the  dispersion  that  tends  to 
spread  out  the  wave  front.   Recently  the  solitary  waves  have  been  discovered 
in  plasma  waves  and  ion-acoustic  waves  as  well  as  in  the  anharmonic 
crystals,  and  much  effort  has  been  concentrated  on  the  analysis  of  such 
waves.   As  for  the  crystal  lattices,  Zabusky  and  Kruskal  showed  that  the 

equation  of  motion  (1.1)  may  be  reduced  to  KdV  equation  or  its  variations 

(9) 
by  introducing  new  coordinates,     and  using  an  appropriate  approximation. 

Toda  found  a  particular  solution  of  the  solitary  wave  type  for  the 

so  called  "exponential  lattice",  which  is  described  by  a  force  function 
f(x)  =  (l-exp(-bx) ) /b .   This  is  the  only  exact  solution  known  for  any  an- 
harmonic lattice. 


As  for  the  ergodic  problems,  many  numerical  experiments  have  been 
done  under  various  conditions.   The  original  FPU  situation  was  re-examined 
by  Tuck  for  longer  periods  of  time,  resulting  in  the  "super  recurrence". 

On  the  other  hand,  the  system  approached  to  the  thermal  equilibrium  for  a 

(13) 
strong  nonlinearity     or  for  an  initial  condition  given  to  the  higher 

(14) 
mode.       The  phenomenon  of  FPU  recurrence  was  studied  from  the  viewpoint 

of  the  normal  mode  coupling  by  Ford,    }'^  Jackson^  '    and  Pasta  e_t 

( 18 ) 
al .      Jackson  could  quantitatively  show  the  behavior  of  a  few  lower 

modes,  which  is  quite  in  good  agreement  with  the  result  obtained  by  FPU. 

The  development  of  the  area  of  the  heat  conduction  will  be  dis- 
cussed more  in  detail  in  the  next  section. 

1.2   Problems  of  Heat  Conduction 

Among  several  published  reports  on  numerical  experiments  of  the 
heat  conduction  problem,  the  one  by  Jackson,  Pasta  and  Waters 

(to  be  referred  to  as  JPW  throughout  this  thesis)  is  very  well  documented 
on  the  computational  aspects  of  their  model  too.   Therefore  we  will  de- 
scribe the  problem  along  the  lines  of  their  report. 

Let  us  consider  a  one  dimensional  lattice  which  consists  of  N 
atoms  of  identical  mass  m  connected  by  anharmonic  springs.   We  choose  free 
boundary  conditions  on  both  ends.   The  force  function  is  of  polynomial 
type: 

f(x)  =  x  -  Xx2  +  "x3  (1.3) 

where  A  and  H  are  the  nonlinear  coefficients  for  the  quadratic  and  cubic 

2  k1 

terms,  respectively.   The  special  relation  h  =  2  \  /3  or  h   =  —     is  used 

to  assure  that  the  force  is  always  attractive.   X  may  be  varied  from  one 

spring  to  another  (JPW  called  this  "defects"),  but  we  will  only  be  concerned 


with  a  constant  \  in  all  cases.   Next  we  define  a  heat  reservoir  in  the 
I'ollowing  way;  every  time  an  end  particle  hits  the  reservoir,  it  loses  its 
velocity  and  is  kicked  back  toward  the  lattice  with  a  new  velocity  v 
sampled  from  the  probability  density  function 

p(v)  =  (v/T)exp(-v2/2T).  (1.4) 

This  definition  was  given  by  JPW  and  is  very  well  suited  for  numerical 
experiments  on  digital  computers.   The  lattice  is  located  between  two  heat 
reservoirs  of  given  temperatures  and  is  subject  to  random  motion.   If  two 
temperatures  are  different  we  expect  an  energy  flow  from  the  high  tempera- 
ture side  to  the  low  temperature  side.   In  other  words  we  can  define  the 
heat  flux  J  numerically  by  the  time-avaeraged  energy  transferred  through 
the  lattice.   We  can  also  define  the  temperature  of  each  atom,  either  by 
time-averaged  kinetic  energy  or  by  time-averaged  total  energy.   Usually  the 
former  is  chosen  for  its  simplicity.   Fourier's  law  of  heat  conduction 
states  that  the  heat  flux  J  is  proportional  to  local  temperature  gradient 
so  that  the  thermal  conductivity  can  be  defined  as  follows: 

J  =-hVT. 

Computationally,  ■"  can  be  calculated  from  the  averaged  heat  flux  J  and  the 
linear-interpolated  temperature  distribution.   JPW  performed  several  types 
of  numerical  experiments  on  the  lattice-reservoir  system  just  described. 
As  they  mentioned  in  the  general  description  of  their  results,  pulse-like 
waves  are  observed  to  propagate  through  the  lattice  without  much  distortion 
in  spite  of  the  irregular  motions  of  atoms  and  defects  (Fig.  2  and  Fig.  3 
of  (21)).   This  fact  suggests  that  the  phenomenon  of  heat  conduction  might 
be  considered  in  terms  of  localized  waves  (such  as  those  we  introduced  in 


jfl.l).   This  leads  to  such  questions  as: 

(1)  Is  it  possible  to  create  solitary  waves  by  impulse  type  initial 
conditions  ? 

(2)  If  so,  do  such  waves  propagate  through  the  lattice  without  appre- 
ciable dispersion  of  energy? 

(3)  How  do  such  waves  contribute  to  energy  transfer? 

(4)  What  are  the  effects  of  boundaries? 

(5)  What  will  be  the  possible  sources  of  the  temperature  gradient  for 
anharmonic  lattices  which  was  observed  by  JPW  and  other,  authors? 

And  so  on. 

Our  main  purpose  in  this  thesis  is  to  establish  a  dynamical 
picture  of  the  heat  conduction  through  the  lattice  in  terms  of  the  localized 
waves  and  also  to  find  the  limitation  ot  such  approach.   We  will  also  dis- 
cuss some  of  the  related  topics  herein. 

1.3   Outline  of  This  Thesis 

This  thesis  consists  of  two  parts,  one  being  on  the  mathematical 
derivations  and  the  other  being  on  the  numerical  experiments.   In  Chapter  2 
we  will  study  the  dynamical  properties  of  the  harmonic  lattices.   An 
emphasis  will  be  put  on  the  solutions  discovered  by  Schrodinger,  which  are 
very  convenient  for  describing  localized  solutions.   Some  of  the  thermal 
properties  of  the  harmonic  lattice  will  be  also  described.   In  Chapter  3 
we  will  discuss  the  properties  of  the  anharmonic  lattices;  first  in  con- 
tinuum approximation  and  then  about  Toda ' s  exponential  lattice.   The  soli- 
tary wave  solutions  and  the  periodic  wave  solutions  are  obtained  in  both 
cases.   The  special  features  of  such  solutions  will  be  extensively  discussed 
and  they  will  be  simply  explained  in  terms  of  elementary  functions.   Chapter 
l\    is  about  the  numerical  methods  and  some  related  problems  encountered  in 


the  actual  computation.   We  will  show  that  Nystrom's  integration  formula 
for  the  special  equations  of  the  second  order  is  quite  advantageous  for  our 
purpose.   In  Chapter  5  we  will  analyse  the  computational  results  on  the 
basic  dynamical  properties  of  our  model  lattice.   Specifically  we  will 
show  that  the  solitary  waves  can  be  created  when  the  lattice  interacts 
with  the  heat  reservoirs.   Since  no  exact  solution  is  known  for  our  model 
lattice,  the  results  obtained  in  Chapter  3  are  applied  whenever  applicable. 
In  Chapter  6  we  will  discuss  the  heat  conduction  problem,  considering  the 
results  in  Chapter  5,  our  own  computational  results  in  the  thermal  situa- 
tion as  well  as  those  obtained  by  JPW  and  other  authors.   Some  of  the 
questions  raised  in  21.2  will  be  answered. 


CHAPTER  2.   THE  HARMONIC  CRYSTAL  LATTICES 

We  consider  the  one  dimensional  crystal  lattices  in  this  chapter 
which  consist  of  either  an  infinite  or  finite  number  of  atoms  of  identical 
mass  and  connecting  springs  which  obey  Hook's  law,  f(x)  =  Kx.   The  equation 
of  motion  lor  the  displacement  y   of  the  nth  atom  from  its  equilibrium 
position  is  then 


d2y 


m  

dt 


Boundary  conditions  and  initial  conditions  are  to  be  described  later. 

4 

The  linear  theory  of  lattice  vibration  has  a  very  long  history 
in  the  theoretical  physics,  or  we  should  rather  say  that  a  considerable 
part  of  mathematical  physics  was  developed  around  this  very  subject  in  its 
early  stages.   Newton  was  the  first  to  investigate  this  problem.   He 
speculated  that  "the  phenomena  of  Nature ••• -may  all  depend  on  certain 
forces  by  which  the  particles  of  bodies* • -are  either  mutually  impelled 
towards  one  another,  and  cohere  in  regular  figures  or  are  repelled..." 
(Principia).   it  is  along  this  line  that  Newton  considered  how  the  sound 
wave  propagates  in  the  air.   By  discretizing  the  continuous  media  such  as 
air  into  particles  of  mass  m  he  obtained  the  equation  of  motion  identical 
to  (2.1).   Note  that  the  theory  of  the  partial  differential  equations  was 
not  developed  at  this  time.   Newton's  idea  was  developed  by  Bernoulli, 
Lagrange,  Cauchy,  Hamilton  and  many  other  well  known  mathematical-physicists. 

For  more  details  on  the  historical  development  of  the  theory,  we  refer  the 

(22 ) 
reader  to  L.  Brillouin.      The  basic  idea  behind  the  development  of  the 

theory  oi  harmonic  lattices  is  the  decomposition  of  the  vibratory  motion 

into  the  Fourier  components,  or  the  normal  modes  of  the  system. 
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2.1   Normal  Mode  Representation  of  the  Vibration  of  the  Harmonic  Lattices 
Suppose  that  we  have  a  finite  number  of  atoms.   The  normal  mode 
representation  is  obtained  by  putting  the  solution  of  (2.1)  in  the  follow 
i  ng  form: 

y  (t)  =  f   exp(iuut). 

Substituting  (2.2)  into  (2.1)  yields  a  finite  difference  equation  for  f  : 

n 

f  J.1-C2-  T^2)f  +f  i  =  °  > 
n+1      K     n   n-1 

/K 

which  we  can  solve  bv  putting  f   =  exp(ikn).    By  defining  U)  =  2*/  —  we 

n  o       jij 

obtain  the  dispersion  relation 

x   =  uu    sin  (k/2)      0  •_'■   k  <   n  .  (2.2) 


The  system  is  dispersive  for  large  k.   The  values  of  possible  k's  depend 

on  the  boundary  condition,  but  in  general  there  are  N  distinct  k's  (we 

call  them  the  modes  of  the  system)  if  there  are  N  atoms  in  the  system. 

The  larger  the  number  of  atoms  in  the  system,  the  greater  the  density  of 

the  allowed  values  along  the  k  axis.   Note  that  uu  is  always  less  than  JJ  . 

o 

This  is  why  the  mass-spring  system  is  often  called  the  mechanical  low-pass 
filter.  Actually  we  can  show  that  there  is  a  close  correspondence  between 
the  mechanical  filters  and  the  electrical  filters  (section  2.3).   Also 


refer  to  L.  Brillouin. 


(22) 


Going  back  to  the  problem  of  finding  the  k's  for  a  given  boundary 

(23 ) 
condition,  we  take  the  following  steps     : 

(i)   Put  f   =  A  exp(ikn)  +  B  exp(-ikn)  and  find  k's  that  satisfy 
the  boundary  condition, 
(ii)   Express  B  (or  A)  in  terms  of  A  (or  B)  to  obtain  a  unique  repre- 
sentation for  f  .   This  gives  a  set  of  orthogonal  functions  from 

n         °  ° 
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which  we  can  obtain  a  set  of  orthonormal  functions, 
(iii)   Obtain  w  for  each  k.   They  are  all  located  on  the  sine  curve 

given  in  (2.2). 

Some  typical  boundary  conditions  are  listed  in  Table  2.1 
together  with  the  values  of  k  and  uu,  and  normal  coordinates. 

There  is  another  type  of  solution  for  an  infinite  lattice  which 
is  obtained  by  putting 

yn(t)  =  fh-exp(-2(3t)  (2.3) 

where 

fn  -  exp(2an).  (2.4) 

Following  the  same  procedures  as  before,  we  can  obtain  the  relation 
between  a.  and  j3  as: 

£3  =  U)   sinh  (a)  .  (2.5) 

Because  of  its  unboundedness  this  solution  is  seldom  mentioned.   If  the 

amplitude  growth  is  suppressed  by  some  nonlinear  effect,  however,  the 

resulting  solution  will  have  an  exponential  behavior  in  the  region  in  which 

the  amplitude  is  sufficiently  small.   This  is  actually  the  case  for  Toda's 

solitary  wave  solution.   We  will  come  back  to  this  problem  in  §3.5. 

Going  back  to  the  normal  mode  representation,  we  can  find  the 

solution  of  (2.1)  for  a  given  initial  condition  by  superposing  the  normal 

modes.   As  an  example,  let  us  consider  the  wave  motion  when  the  left  end 

atom  of  a  free-ended  lattice  is  given  an  impulse  of  unit  magnitude  at  t=0. 

The  general  solution  is: 

p.  B 

y    (t)    =  A  +B    t   +  N  [A   cos(w   t)+  —cos  (a>   t )  ]cos[mn(2n+l) /2N]      (2.6) 

n                  oo  m             mOJ               m 

<-- >  m 
m=l 


Table  2.1 
Normal  Mode  Representation 


periodic 
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y  (t)  =  a  (t)exp(2nimn/N) 


where  a  (t)  =  exp(id)  t) 
m  m 

2    K    2 
and   a)    =  —sin  (TTm/2N) 
m    m 


m  =  1,2,  ■  •  •  ,N 


free  ends 


l(t)  =  a  (t)cos[Tim(2n+l)/2NJ      m  -  051,---,N-1 


where  a  (t)  =  exp(iJj  t) 
m  m  ' 

2        K        2 
and        a,        =  —sin    (TTm/2N) 
m  m 


Oy*5K>TnfK>wK>wxiD 


fixed    ends 


V    (t)    =   a    (t)sin(rimn/N) 
'  n  m 


m    =    1,2, • • • ,N-1 


where   a    (t)  =   exp(iO-    t) 
m  r         m 

9  K        2 

and        <k        =  -sin    ("m/2N) 

m  m 


^T<>^!I^O^^<>^^>yTVt 
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where  A  and  B  are  the  Fourier  coefficients  of  the  initial  displacement 

and  velocity  which  will  be  given  shortly.   We  can  also  obtain  y  (t)  as 

^n 

follows : 

N-l 

Yn(t)    =   Bq+  V)       [Bmcos(uJmt)-AmiUmsin(u)mt)]cos[mn(2n+l)/2N].  (2.7) 

m=l 


Initial  conditions  are 


N-l 
\ 


y  (0)  =  A  +   ;■   A  cos[mrr(2n+l)/2N]  (2.8) 

n       o       m  v    ' 


m=l 

3nd  N-l 


y  (0)  =  B  +  ;   B  cos[mH(2n+l) /2N] .  (2.9) 

n       o       m 


m=l 


Or  solving  for  A   and  B  , 
m      m' 


N-l 

e   ,  -, 

A   =  -7  S   y  (0)cos[mn(2n+l)/2N]  (2.10) 

m    N       n 


n=0 

B      = 
m 

N-l 
e 
m 

N 

y  (0)cos[mn(2n+l)/2N]  (2.11) 

n=0 

,         . 1  i  f  m=0    .     ,  _  . 

where  E   =-   .      .   A  and  B  give  the  motion  of  the  center  of  mass  of 
m    2  if  m?t0     o       o  ° 

the  lattice.   For  the  present  problem,  we  have  the  following  initial 
condi  tions : 

y   =0  for  all  n   —  -rr—  ;■  A   =  0  for  all  m, 
n  m 

E 

y   =  6  n  >  B  =  -7cos(mn/2N)  . 

^n     n,0  /      m         N 

Substituting  A  and  B   into  (2.7)  yields 

m  m 

N-l 
y    (t)    =    1/N  +   2/N  cos(mn/2N)cosLtsin(mn/2N)]cos[m(2n+ln/2N)].     (2.12) 

m=l 
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This  is  the  solution  for  the  given  initial  condition.   Note  that  the 

amplitude  distribution  among  the  normal  modes  is  not  uniform  but  shows  a 

monotonous  decrease  as  the  mode  number  increases.   The  energy  distribution 

is  obtained  from  the  square  of  the  amplitude  distribution,  since  initially 

there  is  no  potential  energy  and  the  total  energy  of  each  mode  is  strictly 

conserved.   Both  distributions  are  shown  in  Fig.  2.1  and  Fig.  2.2.   For 

the  anharmonic  lattices  with  the  same  boundary  and  initial  condition  we 

can  also  use  this  distribution  approximately  if  time  is  very  close  to  t=0. 

This  will  be  discussed  in   §5.1. 

Next  let  us  consider  the  case  when  N  goes  to  infinity.   From 

(2.12)  we  have 

4   >n/2 
v  (t)  =  —      cos(x)cos[t  sin(x)  ]cos[  (2n+l)xJ  dx  +0(1/N) 
-  n 


TT    J 

o 


(N  -»  °°) 


where 


where 


mrr         n 
X  =  2^  '  dx  =  2^  ' 

4  ^ 
=   ~  J   cos(wt)  cosL (2n+l )arcsin(w) J  dw 


w  =  sin  x,  dw  =  cosx  dx  , 

1 
=  "|   j  cos(wt)(-l)n  U2n+1(w)  dw,  (2.13) 


where  U    1 (w)  is  the  second  kind  of  solutions  of  Tchebychef f ' s  equation: 

2   2        2  2 

(1-w  )d  U  (w)/dw  -wdU  /dw+n  U  (w)  -   0.   Making  use  of  the  following  identity 
n           n       n 

for  an  odd  n,  .  .  ,  ,  , 

o   1        ,  U  (w)  J  (t) 

-   j   (-i)  ~ cos(wt)  dw  =  — - —  ,         (A. 7) 

o 


See  Appendix  A  for  the  derivation. 


AmpMuJe 


l"-- 
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Fig.  2.1.   Amplitude  distribution  (N=9) 
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Fig.  2.2.   Energy  distribution  (N=9) 
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we  can  obtain  the  expression  for  y  (t)  when  N 


J2n+l(t) 


y  (t)  =  2(2n+l) .  (2.14) 

n  t 


We  will  obtain  this  solution  much  more  easily  in  the  next  section.   This 
solution  is  the  response  of  the  velocity  of  the  nttri  atom  to  a  unit  impulse 
at  the  left  end.;  it  shows  how  the  energy  propagates  through  the  harmonic 
lattice.   In  a  thermal  situation  initial  impulses  are  given  at  random. 
We  will  come  back  to  this  problem  in  §2.4. 

2.2   Schrodinger ' s  Solutions 

A  different  type  of  solutions  was  found  by  Schrodinger  in  1914 

(24) 
for  infinite  harmonic  lattices.       Since  his  solutions  are  very  convenient 

for  describing  the  energy  transport  through  the  lattice  we  will  summarize 

some  of  their  features. 

Let  us  introduce  the  new  coordinates  Tj   according  to  Schrodinger; 


TL   -  v'm  y  (2.15a) 

2n       n 

\m  -  ^<vW  .  (2-:5b) 

2 
and  t'  =  u;  t,  where  uo    =  4K/m.    Then  (2.1)  yields 
o  o 

dTL  /dt'  =  (\        -7)    )/2  (2.16a) 

2n         2n-l   2n+l 

dT     /dt'  =  (T|   -T]     )/2  .  (2.16b) 

2n+l  2n   2n+2 


These  equations  are  equivalent  to 


dT,  /dt'  =  (Tl   -71  ,)/2  (2.17) 

n         n-1   n+1 

u-  t.  ■  (25) 

which  is  a  recurrence  formula  for  Bessel  functions  and  Neumann  functions. 
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Since  Bessel  functions  with  integer  order  are  the  only  ones  which  are 
finite  for  all  t'  >  0,  we  will  restrict  ourselves  to  such  solutions. 
Therefore  the  general  solution  is  given  as  a  linear  combination  of  Bessel 
functions  as  follows: 


k 

where  a,  's  are  determined  from  the  initial  conditions.   Let  us  consider 
k 

an  initial  impulse  given  to  the  0    atom  when  all  the  atoms  are  in  their 
equilibrium  positions.   Eventually  we  will  consider  the  example  discussed 
in  §2.1  in  terms  of  Schrodinger ' s  solutions.   We  have,  for  the  present 


example,  a.  =  vm  &,  ^ ,  or 
k       k,0 


T]   =  Vm  J(f) 
n        n 


Going  back  to  the  original  coordinates,  we  have 


y„  -  Jo  ^o^  (2-L9) 

n    2n   o 


Vyn+1    ■    ''WK    J2n+l(V> 

-i"J2»rt<V>-  (2-20) 

O 

In  general  Schrodinger ' s  solutions  are  very  convenient  for  studying  the 
propagation  of  local  disturbances  through  the  lattices.   Fig.  2.3   shows 
how  an  initially  localized  disturbance  propagates  with  time.   As  is  clear 
from  Fig.  2.3  ,  the  harmonic  lattice  is  very  dispersive  and  energy  spreads 

out  rather  quickly.   From  the  asymptotic  formula  for  Bessel  functions,  the 

,    ,"1/2  (25) 
amplitude  of  the  wave  is  known  to  decay  as  (w  t; 

In  order  to  satisfy  the  free  boundary  condition  at  the  left  end 

of  the  lattice,  we  add  another  impulse  of  the  same  magnitude  to  the  atom 
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19 
on  the  left  of  atom  n=0.   Then  the  solution  is 

*n  =  J2n(V}  +  "W^o0  *  (2'21) 

In  particular  the  atom  n=0  and  n=-l  oscillate  exactly  in  the  same  way. 
Initially  there  is  no  energy  stored  in  the  spring  that  connects  them; 
therefore  atom  n=0  actually  feels  no  forces  from  its  left  neighbor.   Con- 
sequently we  can  eliminate  all  the  atoms  with  negative  indices,  and  the 
resulting  lattice  is  semi-infinite.   By  making  use  of  a  recurrence  formula 

<Jn-l+Jn+l)/2  '  "  VX'  <2'22> 

(2.21)  may  be  written  as 


y   =  2(2n+l)  J_  J_.  (cu  t)/w  t  .  (2.23) 

n  zn+1   o    o 

We  can  also  obtain  the  response  of  y  -y   ,  to  a  unit  impulse  in  the  same 

r  J  n      n+1  r 

way : 

y   -y    ,,    =   2/oj   [j_      .  (J    t)    +  J.      -(cu   t)] 
yn   ^n+1  '    o      2n+lv   o  2n+3v   o 

=   8(n+l)/x   [j  (a,   t)/0)  t]    .  (2.24) 

o      2n+2      o  o 

Fig.  2.3  shows  how  each  atom  responds  to  an  impulse;  we  can  see  that  each 
wave  decays  somewhat  faster  than  in  the  previous  example.   This  is  due  to 

the  fact  that  there  is  no  force  to  pull  atom  n=0  back  to  the  original 

-3/2 
position.  (2.24)  asymptotically  behaves  as  (u>  t)     .   In  the  similar  way 

we  can  obtain  the  velocity  response  function  when  the  initial  impulse  is 

applied  to  the  kth  atom.   The  results  are 

*„«>  "  Ja(k-n>(V>  +  J2(n+k)+2(V>       0<n<K 
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yk(t)    =    Jo(a,ot)    +   J4k+20V>  n    =k 

^n(t)    =    ^(n-k)^    +   ^(n+WV*  "    >k    ' 

It  is  easy  to  see  that  the  free  boundary  condition  is  satisfied  at  n=0. 

Parenthetically,  the  solution  we  just  obtained  is  also  a  solution 
to  the  infinitely  long  lattice  with  atom  n=0  twice  as  heavy  as  the  other 
atoms,  to  which  an  impulse  is  applied. 

Schrodinger ' s  solutions  are  very  simple  and  clear  as  long  as  the 
lattice  is  infinite,  but  in  the  case  of  a  finite  lattice  the  reflection  of 
waves  must  be  taken  into  account.   Consequently  we  have  an  infinite  series 
of  Bessel  functions  instead  of  a  single  term,  even  for  the  response  func- 
tion to  an  impulse.   For  example,  the  velocity  impulse  function  for  a 
lattice  with  N  atoms  (n=0  and  n=n-l  being  free  ends)  is: 


F  (oj  t)  -  J0  (U)  t)+J.  ,„0  t)+J       (ou  t)+J.M    (uj  t) 
n   o      2n   o    2n+2   o    4N-2-2n   o    4N-2n   o 


+  ^N^ntV^WWV^SN^^n^o0 


+JQM  0  (a)  t) 

8N-2nN  o 


+ 


=  ^J[j4kN+2n(u,ot)+J4kN+2n+2(Uiot)+J4(k+l)N-2n^ot) 
k=0 

+  J4(k+1    -2<V>]-  (2'25) 

We  can  easily  show  that  the  free  boundary  conditions  are  satisfied  (the 
method  of  mirror  images).   The  response  function  (2.25)  is  identical  to  the 
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solution  which  we  obtained  in  §2.1,  (2.12).   As  we  can  see  from  this  example, 
the  normal  mode  representation  and  Schrodinger ' s  solutions  are  complimentary 
to  each  other. 

2.3  An  Analogy  between  the  Harmonic  Lattices  and  Electrical  Circuits 

There  is  a  one-to-one  correspondence  between  the  harmonic 

lattice  and  electrical  low  pass  filter  (L-C  ladder  circuit)  as  has  been 

(22) 
discussed  by  many  authors,  for  example  L.  Brillouin.       It  is  often  very 

convenient  to  consider  an  equivalent  electrical  circuit  of  a  dynamical 
system;  especially  the  notion  of  impedance  in  the  former  enables  us  to 
skip  some  of  the  lengthy  calculation.   As  an  example  we  will  consider 
Schrodinger ' s  solution  here.   Table  2.2  summarizes  the  correspondence. 
The  normal  modes  of  an  electric  circuit  such  as  a  low  pass  filter  may  be 
obtained  by  looking  at  the  zeroes  of  the  impedance.   Hence  the  normal  modes 
of  the  corresponding  dynamical  system  are  easily  obtained  by  simple  substi- 
tution.  The  impedance  is  also  very  useful  when  we  consider  the  wave 
reflection  at  an  impurity  in  the  harmonic  lattices.   This  is  a  problem  of 

impedance  matching  in  electrical  circuits.   But  we  will  not  go  any  further 

(  7  ft  ^ 
in  detail  on  this  subject.   Refer  to  van  der  Pal 

2.4  Heat  Flux  and  Temperature  Calculation  for  Semi-infinite  Harmonic 
Lattices 

As  a  simple  application  of  Schrodinger ' s  solutions  we  calculate 
the  heat  flux  and  the  temperature  distribution  of  the  semi-infinite  har- 
monic lattices  which  we  discussed  in  §2.2.   We  suppose  that  only  the  atom 
at  the  end  of  the  lattice  interacts  with  the  heat  reservoir;  it  receives 

random  impulses  from  the  reservoir.   The  approach  we  take  here  follows 

(27) 
ol  niilenbeck  and  Ornstein  in  their  classical  work  on  Brownian  motion. 


Table  2.2 
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Dynamical  System 


Electric  Circuit 


difference  of  the  displace-    \    V  :   voltage  of  the  n    capacitor 
merit  of  two  neighboring 
atoms  (i.e. ,  yn_1Yn) 


S   = 


1  dyn  ' 

:   velocity  of  n*-"  atom    I 

K   dt  J  ;   n 


current  through  nth  inductance 


m:    mass  of  the  atom 


L:    inductance 


k:    spring  constant 

boundary 

free  end  on  the  left 


C   :  (reciprocal  of)  capacitance 


short  circuit  on  the  left 


a)  =  2v'K/m 
o 


1  0)  =  2/v/LC 
o 


initial  condition 


yo(t)  =   VQ/2K6(t) 


!  V  =  E  6(t) 
(   o    o 


equations 


1/Kdu  /dt  =  S   . -S 
n       n-1   n 


m  dS  /dt  =  u  -u   n 
n       n  n+1 


solution 


S  (t)  =  v  /KJ_  (uj  t) 
n       o    2n   o 


CdV  /dt  =1    -I 
!     n       n-1   n 


Ldl  /dt  =  V  -V  ,n 
n       n  n+1 


I  (t) 
I   n 


2E  /L  J_  (oj  t) 
o     2n   o 


or 


yn^t}  =  VoJ2n<V> 


Table  2.2  (continued) 
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impedance 


.Uj 


Z(io)    =  vmK  J     L-(— ) 

o 


Z(uu)    = 


c 


?o         <J|  ^         !fe  ^ 


m 


w       m       m       m 


K       K       K        K 


Jg.     _*i        i*        x3 


^4  Z, 


*n         c      c      c      c 

-*JL         7J77       tJtt       -Xr       -Xt 
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We  have  obtained  the  velocity  response  function  in  §2.2: 


F  (a)  t)  =  J   0»  t)+J     (ou  t) 

no       2n  o     Zn+Z   o 


=  2(2n+l)  J2n+I(^0t)/^ot.  (2.23) 


The  velocity  of  the  n   atom  at  time  t  is  therefore  given  as 


(t)  =  P   f(T)  F  (to  (t-T))dT  (2.26) 


y 
n      °  no 


where  f(T)  is  the  random  function  which  is  characterized  by  f(T)=0   and 


f(x)  f(r')  =  A6(t-t').   (    )  means  the  phase  average  of  the  variable 
given  in  the  parentheses.   Since  the  effect  of  any  initial  condition  is 
supposed  to  have  propagated  toward  n=+°°  by  the  time  we  observe  such 

ensemble  of  lattices  at  time  t  (actually  amplitude  and  velocity  of  each 

-1/2 
atom  damp  approximately  as  (uu  t)    ),  we  only  consider  the  response  of  each 

atom  of  lattices  due  to  the  random  forces  acting  on  the  0   atom  during 

-2  2 
the  time  interval  (-co,  t)  .   A  constant  A  has  the  dimension  of  (LT   )   and 

is  to  be  determined  later  to  make  the  temperature  of  the  end  atom  kT/2m. 
By  using  an  explicit  form  of  F  (t)  in  (2.23),  we  can  calculate  the  phase- 
averaged  squared-veloci ty  for  each  n  as: 

r  t  t 

y    (t)      =  dT  dT1    f(T)f(r')    F    (uu    (t-T))F    (uu    (t-T1)) 

n  J  "  no  no 

_  00  -00 

9  ? 

=    A/UU  F     (CD     (t-T))d(JL)     T 

O         J  O  O 

-  00 

CO 

2     -  2 

=  A/cu         I      F      (s)    ds  (s   =   UJ    (t-T)) 

o      "J         n  o 

4/      2   4°2n+1>2   JL+1     , 

=  A/uu ds 

o  z 


25 


4A(2n+l)2    r(2n+l/2) 
(^o222r(3/2)r(3/2)r(2n+l+3/2) 

A(2n+l)2/ra   2 
o 


=  2 ,  (2.27) 

(2n+l)    -1/4 

where  the  phase  average  is  taken  at  time  t  over  an  ensemble  of  lattices 
(a  large  number  of  similar  but  independent  lattices),  which  have  started 

at  t=-oc  with  identical  initial  conditions  given  to  the  corresponding  atoms 

(25) 
therein,  and  we  have  also  made  use  of  the  following  formula 


J   Jm(°  V°     m   Y(\)    ((m+n-W/2) 


°      t  2K    r((\4tn-n+l)/2)  r((\+n-m+l)/2)  f(  (m+n+X+1)  /2) 


(2.28) 


We  can  determine  the  constant  A  from  the  condition  y  (t)  =  kT/m,  that  is, 

2 
A  =  3hju  kT/4m   or 
o 


f(T)f(T')  =  3thju  kT/4m6(i-T')  . 
o 

Fig.  2.4  shows  the  temperature  distribution.   There  is  a  slight 
decrease  of  the  temperature  near  the  end  atom,  then  the  temperature  becomes 
almost  uniform.   This  tendency  is  not  changed  if  we  define  the  temperature 
of  each  atom  by  the  total  energy  instead  of  kinetic  energy.   The  heat  flux 

(instantaneous  flux  at  each  point  of  the  lattice)  is  obtained  by  using  the 

(21) 
following  formula     ; 


JH(n)  -  -K(yn+yn+1>(yn+1-V/2  ■  (2'29) 


We  can  obtain  the  first  term  immediately  from  (2.23): 

(VW/2  ■  f    f('»J2„("ot)+2J2n+2("ot)+J2„rt(a,ot))dl/2-         (2-30) 
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it.,  m 


n 


Fig.  2.4.   Temperature  distribution  for 
semi-infinite  lattice. 
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To  calculate  the  second  term  we  have  to  define  the  response  function  G  (t) 

n 


for   y    ,,-y    .      From    (2.24)    we   define 
yn+l      n 


G    O   t)    =   -       N        ^        """• — - —      .  (2  3D 

no  U)  a)    t  ^ .  ~  l  ; 

o  o 


((D   t)    =   .    8(n+l)      ^n+Z^o^ 

O  'X1  O)    t 

O  o 

Unlike   F    (u)   t)    this    function   has    the   dimension   of   time.      Using  G    (a   t)    we 
n      °  no 

can   obtain   the    second    term  as   a   convolution   form: 

t 
Vn^o0    "    yn(V}    =  J"      f(T)Gn(^(t-T))dT.  (2.32) 

—  OC 

Substituting  (2.30)  and  (2.32)  into  (2.29)  yields 

mSlaa3*(f  f(T)(J2n(«»ot)+2J2n+2(«.ot)+J2      (.  t))dT) 

o       _ot> 

X  Cl   f(^>     '  (°_T) dT')  (2.33) 

where  we  have  neglected  the  effect  of  initial  conditions  as  before.   Taking 

the  phase  average  and  changing  the  variable  x  (t-T)=s  yield 

,,    .,  v    ...      a     J0    (s)J0   ^(s)  2  J.    ,0(s)        J0   ^0(s)J0   ^  (s) 

~ — 7~"T  4  (n-fl)aAK    p      ,   2n  2n+2  2n+2  2n+2  2n+4 

Ju(n)      =  — ^ 1 ( )    +  +  ds 

H  3  "  s  s  s 

ju  o 

°  2 

_   6nkT(n+l)K     •        2n+2^' 


■  ds 

ma  J  s 

o  o 


3nkTK 

2  ma 
o 


3nkTai 
o 


8 


(2.34) 


where   we    have   made    use   of    the    formula    (2.28)    for    \=l ,    or   more   precisely, 

"  Jm(°    Jn(t)  2    sin((m-n)n/2)         A         .  , 

dt    =  *"r t1 — ' —     =0        if   m-n:      even 

o  (m   -n    )  tt 

and 
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«  J  2(t) 

^     m        At-  —L 

J   ~1     dt  =  2m  • 
o 

Hence  a  constant  heat  flux  exists.   Since  the  temperature  distribution 

becomes  more  and  more  uniform  as  n  increases,  it  is  obvious  that  Fourier's 

law  of  heat  conduction  does  not  hold  in  this  case:   we  cannot  define  the 

thermal  conductivity!   Similar  results  have  been  obtained  by  several  authors 

/TO  \    /O  Q  \ 

for  different  boundary  conditions  and  different  thermal  situations. 

The  results  obtained  in  this  section  entirely  depends  on  the  fact 

CO 

r    2 
that  the  integral  (   F   (t)dt  converges.   That  is,  the  energy  fed  by  the 

o 

random  force  to  the  left  atom  is  transferred  to  other  atoms  at  a  proper 

rate.   As  for  classical  problem  of  the  Brownian  particle  (Langevan  equation) 

the  energy  fed  from  the  random  force  is  dissipated  by  the  viscosity  so 

that  there  is  no  accumulation  of  energy  in  the  Brownian  particle.   If  we 

consider  the  Schrodinger ' s  lattice,  for  which  the  velocity  response  func- 

f*  h 

tion  is  given  as  F  (t)  =  J   (cu  t)  (suppose  that  the  0   atom  interacts  with 
the  heat  reservoir),  the  temperature  of  each  atom  rises  indefinitely  since 

00 

r        2 
the  integral  (   J_  (^  t)dt  diverges  for  all  n.   Such  catastrophe  also 

o 

occurs  for  the  harmonic  lattice  of  finite  length  if  there  is  no  dissipation 

of  energy. 
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CHAPTER  3.   THE  ANHARMONIC  LATTICES --CONTINUUM  APPROXIMATION 
AND  TODA'S  EXACT  SOLUTIONS 


We  can  approximate  the  system  of  nonlinear  ordinary  differential 
equations  (1.1)  by  a  nonlinear  partial  differential  equation,  provided 
that  the  relative  displacement  of  any  two  neighboring  atoms  is  not  large. 
In  the  first  part  of  this  chapter  we  will  treat  such  approximation  and 
summarize  two  special  solutions.   These  solutions  have  their  complete 
analogues  in  the  water  wave  theory ,  in  which  they  are  called  the  cnoidal 
wave  (periodic  solution)    and  the  solitary  wave  (localized  solution).   ' 

A  new  name  "Soliton"  is  frequently  used  for  the  latter  because  of  its 

(4  ) 
particle-like  behavior  (due  to  Zabusky     ). 

We  start  from  the  equation  of  motion 

,2 

y„  =  f(y^,-y J  -  f(y  -y„_,>  •  (3-D 


,  2  Jn  VJn+l  'n'     v/n  'n-1 

dt 

We  assume  that  f(x)  ~  x  when  x  ~  0,  that  is,  when  the  lattice  is  near  its 

natural  length. 

Historically,  Zabusky  and  Kruskal  were  the  first  to  introduce 

(  9  ) 
the  continuum  approximation  to  (3.1).      They  reduced  (3.1)  to  a  partial 

differential  equation  for  y(x,t),  which  is  equivalent  to  (3.5).   Then 
they  reduced  it  into  the  KdV  equation  by  considering  the  waves  which 
propagate  in  one  direction  (say,  in  the  positive  direction).   Since  the 
equation  (3.5)  has  solutions  which  are  similar  to  those  in  the  KdV  equa- 
tion and  (3.5)  can  also  describe  the  waves  which  travel  in  the  opposite 
direction,  we  will  use  (3.5)  as  our  basic  equation.   Recently,  Wadati  and 
Toda  showed  that  the  two  solitary  wave  solutions  of  both  equations  give  a 
very  good  agreement.   Readers  are  referred  to  (30)  and  references  therein 
for  the  KdV  equation. 
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In  §3.1,  §3.2  and  §3.3,  we  will  consider  the  polynomial  type 
force  function,  specifically  the  quadratic  nonlinearity  and  the  cubic  non- 
linearity,  for  which  we  can  find  the  analytic  solutions.   We  will  show  that 
both  the  periodic  wave  and  the  solitary  wave  can  exist.   For  the  heat  con- 
duction problem,  the  periodic  solution  is  quite  irrelevant.   But  it  will 
also  be  treated  in  full  since  it  seems  to  be  the  fundamental  solution,  of 
which  the  solitary  wave  is  a  special  case.   For  example,  only  the  periodic 
wave  can  exist  as  a  travelling  wave  solution  under  certain  conditions. 

In  §3.4  we  will  simplify  the  force  function  by  introducing  the 
piecewise-linear  approximation.   We  will  show  that  the  essential  features 
of  the  solutions  obtained  in  §3.2  and  §3.3  are  still  retained  in  this 
approximation.   By  this  model  we  can  explain  why  and  when  the  solitary  wave 
solution  is  possible. 

The  contents  of  §3.4  has  not  been  discussed  before.   In  §3.5 
we  will  summarize  Toda ' s  solutions  for  the  exponential  lattice. 


3.1   Polynomial  Type  Force  Function 

In  the  following  sections  we  will  consider  a  cubic  polynomial 

f  (x)  =  x  +  ax2  +  £x3  (3.2) 

for  the  force  function.   It  is  more  convenient  to  introduce  the  relative 

displacement  of  two  neighboring  atoms  u  =  v   ,  -y   instead  of  y   itself. 

°        n    -  n+1  J n  n 


Then  the  equation  of  motion  for  u   is  obtained  from  (3.1)  as 

H2 

-^u  =  f(u  ..)  -  2f(u  )  +  f(u    )  .  (3.3) 

,  I      n       n+1         n       n-1 
at 

We  will  denote  the  spatial  coordinate  by  x  in  the  continuum  approximation. 
Then  we  have 

,u/du\     h  i°   un    ,  "  /o  uN     n  ,ou, 

Vi  =  ux=n±h <a?     +  2-(T2>     ±  T^T$     +  2T(74}     +'  •  • 

—  x=n  ox      x=n  ox      x=n  ox      x=n  ,_    , N 

(3.4) 
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where  h  is  the  interatomic  distance  in  equilibrium.   Substituting  (3.4) 
into  (3.3)  with  (3.2)  yields 

u  =  h  u"  +  2oh2(uu1,+u*2)  +  3ph2(u2u"+2uu*2) 
.4 
+  -u«"+  ••• 

where  •  means  the  time  derivative  and  '  means  the  spatial  derivative. 

Putting  h2  =  C  2,  t :C  2  =  2ah2,  36h2  =  p-C  2  and  h  =  £-  =  -2-  , 
o     o  ^       o  12   12   ' 

u  =  C  2[u"  +  e(uu')'  +  u.(u2u')']  +  hu,,m  .        (3.5) 

In  order  to  consider  the  travelling  wave  solution,  we  put 

u(x,t)  =  w(x-vt)  -  w(z) 

where  v  is  the  propagation  velocity  of  the  wave  (the  phase  velocity  in 
case  of  periodic  solution)  which  is  to  be  determined  in  relation  to  the 

amplitude  and  other  factors.   Then  (3.5)  yields 

2  2 

?        ?  *C  u.C 

(v  -C   )w"  =  hw""  +  — 2_  (wz)m  +  — 2_  (wj)m  _ 
o  L  J 

Here  '  means  the  derivative  with  respect  to  z.   Integrating  both  hand  sides 

twice  yields 

2         2       2    2 

•  ■       o    3      o    2  ,      o      C     D' 

w"  =  -  -r w   -  — w   +  w  -  —  z  -  —       (3.6) 

3h        2 'i  ^       K     K       v 

where  C  and  D'  are  the  integration  constants.   We  eliminate  the  term 

C 
-—  z  since  it  contains  z  and  is  regarded  as  a  spurious  term  which  appeared 

in  the  process  of  integration.   This  can  also  be  verified  by  approximately 

(3.1)  by  a  nonlinear  partial  differential  equation  for  y  and  comparing  its 

integrated  form  (in  which  w  corresponds  to  -tp(z))  with  (3.6).   So  we  adopt 

(3.6)  with  C  =  0.   Multiplying  (3.6)  by  w1  and  integrating  with  reg.ard  to 
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z  once  yield 

2  2        2    2 

9      hlC         EC        v  -C     . 
,  ,  N2        o    4      o    J        o    2    „ 

(w  }  =  "  ~eZ~  w  "  T5T  w  + — ~w  +  Dw  +  E         (3-7) 

2D' 
where  D  =  -  and  E  =  2E'  are  the  new  integration  constants.   The  right 

hand  side  of  (3.7)  will  be  denoted  by  P(w).   We  will  find  the  solutions  of 

(3.7)  in  the  following  sections.   An  emphasis  will  be  put,  whenever  the 

periodic  solution  is  obtained,  on  the  wave  for  which  the  average  over  its 

period  vanishes  (w  =  0) .   Such  wave  motion  takes  place  around  the  equilibrium 

positions  of  atoms  when  the  crystal  lattice  is  of  natural  length.   If  the 

external  force  is  applied  to  the  lattice  so  that  it  is  either  compressed  or 

expanded,  a  periodic  wave  is  set  on  around  the  new  equilibrium  length.   This 

corresponds  to  considering  a  solution  for  which  w  ^  0;  w  >  0  for  the  expanded 

lattice  (to  be  called  "positive  background")  and  w  <  0  for  the  compressed 

lattice  ("negative  background") .   The  propagation  velocity  of  the  waves 

under  such  circumstances  differs  from  that  of  the  waves  in  the  lattice  of 

natural  length  due  to  the  nonlinearity  of  the  force  function. 

As  for  the  solitary  wave  solution,  we  are  most  interested  in  the 

one  for  which  the  lattice  is  of  natural  length  at  z  =  +  °°.   This  means 

that  the  energy  is  concentrated  at  the  wave  front.   (Fig.  3.1)   We  call 

such  wave  zero-back  ground  wave. 


3.2   Quadratic  Nonlinearity  (u.=0,  e^O) 

The  equation  (3.8)  yields,  in  this  case 

2       2    2 

dw  2      KCo    3    V    o    2 
(^)   =  "  -3^-  w  +  -f-   w  +  Dw  +  E  =   P(w).  (3.8) 


33 


equi'libn'm 


I 


\ x^  \  \  \  \  \  \  ' 

\  \  \  v  \  N  v  V 

\  \  \  \  \  \  \  \ 


— ^       30  h  tatv  u/q.  ve. 


Fig.  3.1.   Solitary  wave  on  zero  background, 


34 


3.2.1   e  >  0 

Periodic  Wave 

We  consider  the  case  in  which  P(w)  has  three  distinct  real  roots: 

2 
P(w)  =  P  (a-w) (w-b) (w-c)     (a  >  b  >  c) 

2  2    2 

r,  EC   ^      _  V   "C  L 

I             o    2  ,   ,   N        o 
p   =— ,P  (a+b+c) —   . 

Then  we  have  a  bounded  solution  in  [b,a]  (Fig.  3.2)  which  may  be 
expressed  in  terms  of  the  Jacobian  elliptic  functions  as 


where 


w(z)  =  a-(a-b)sn  \~~^ —  Pz,k)  =  c+(a-c)dn  (^— ^ —  pz,k) 

9     o  —  K  •>'•■ 

where  k   =  is  the  modulus  of  the  elliptic  function.     Going  back  to 

a-c  r  ° 

x,t  coordinates,  we  have 

2  2K 
u(x-vt)  a  a-Asn  ( — (k  x-u)t),k)  (3.9a) 

=  c+  ^   dn2(^(kxx-cut),k)  (3.9b) 

k 

where 


and 


A  =  a-b 

(amplitude 

TTC 

'  3h 

kx  ""  4Kk  v 

(wave  number) 

u)  =  vkx 

(frequency) 

K  =  K(k) 

(the  complete 

first  kind)  . 


The  velocity  v  is  given  in  terms  of  a,b,c  from  the  quadratic  term  in  (3.8): 

2    2  2 

V   o     2  o 

jf-  =    P  (a+b+c)  =  —^~   (a+b+c) 

or 

v2  =  Co2(l  +  j  (a+b+c))  .  (3.10) 


See  Appendix  B  for  collections  of  solutions 
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.  3.2.   Quadratic  nonlineari ty  (e  >  0) :   periodic  wave 
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Note  that  v  >  C   if  a+b+c  >  0  and  vice  versa.   The  periodic  solution  is 

o 

completely  determined  by  three  zeroes  of  P(w).   If  we  further  assume  the 
zero  background  wave,  we  have  the  two  parameter  family  solutions,  for  which 
A,  k  ,  uo  and  k  are  explicitly  related.   First  we  should  note  that 


-^  J   sn2xdx  =  -^  J   -|(l-dn2x)dx  =  -|(1  -  |^) 


2K 


2K  J    .  2V      y     .  2 
ok  k 


K(k) 


n  2  E  (k ) 

since  Z(x)  =  Jdxdn  x  -   ^  x  is  periodic  of  period  2K  (Zeta  elliptic  func- 

o 

tion  of  Jacobi) .   Here  E(k)  is  the  complete  elliptic  integral  of  the  second 


kind.   Using  the  above  relation,  we  have 


or 


Therefore 


u  =  a-A  sn  =  a-  — r(l-  — )  =  0 
k2    K 


A  n      E^ 
a  =  -2<1-  -) 


u(x-vt)  =  — (dn  (—  (kx  x-u)t))  -  -) 
k 


A_  d_ 
,  2  dZ 


Z(z) 


z=  —  (k   x-a)t) 
rr  x 


From  the  definition  of  A  and  k  ,  we  also  have 


b  =  a-A 


c  =  a- 


2  " 


Substituting  a,b,c  into  (3.10)  yields 


2     ?  e  A 

v  =  C  Z(l+  ^(3a-A-  ^r)) 

k 

=  c  2(1+1(3A  _  3AE  _ 

o  U+  3\2   ,2V   A 

k     k  K 

r1(^+   £A  a   i2   3E 
=  Co  (1+  „2(2_k  "   K)} 


k2>) 


(3.11) 
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Hence  the  velocity  of  the  periodic  waves  can  become  either  "supersonic" 

(v  >  C  )  or  "subsonic"  (v  <  C  )  depending  on  the  degree  of  deformation 
o  o 

(i.e.  the  modulus  k) .   The  critical  modulus  k   is  given  by  solving 

3E(k)      2 
■  ' ),  (  =  2 - k  ,  for  which  v  =  C   holds  independent  of  the  amplitude  A.   Approxi' 
K(k)  o  r  r  ^ 


mately  k  ^  .  98. 
J      c 

The   dispersion   relation    for    zero   background   wave    is 

2  2,     2        „   2  ,,       £    N1     2  AK   2    ,3E      _N1     4 

uu=vk        =  C      (1-   -A)k        -    k(— -)       (— -  -2)k 
x  o  3        x  n/vK  x 

A  2     2   2     4 

If  k  »  0  and  — r  finite,  'Jj   =  C   k   -4nk    (dispersive  harmonic  wave)  and 
2  o   x      x      r 

if  k  >  1  uu   =  C   k   (H — r)  (approaches  to  solitary  wave). 

o   x     3 

Solitary  Wave 

If  b  approaches  to  c,  we  have  the  solitary  wave  solution.   Since 

k   =  •  1 ,  dn   +  sech  .   In  this  limit  (3.9b)  reduces  to 

a-c 

2 
u(x-vt)  =  b  +A  sech  (ctx-pt)  (3.12a) 

where  A  =  a-b 


v'a  -  b   .   /  >--A 
a  =  J 


3*    V  12h 

v2  =  Cq2(1+  ^(a+2b))  (3.12b) 

and  3  =  u*v  . 

Consider  the  zero-background  solitary  wave,  i.e.  b  =  0.   Then 

2 
u(x-vt)  =  Asech  (ax-pt)  (3.13a) 


where 


v2  =  C  2(1+  |A)  (3.13b) 

o     J 

and  \i   =  a*  v  . 

Note  that  the  solitary  wave  is  always  created  toward  the  positive  w  direc- 

i  Lon.   Since  u   =  y   , "V   >  0,  this  means  that  only  rarefactive  solitary 
n    ^n+1  Jn 

wave  is  possible  for  t    >   0. 
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Fig.  3.3.   Quadratic  nonlinearity  (£  >  0)  :   solitary  wave 


Fig.  3.4.   A  typical  profile  of  solitary  wave. 
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3.2.2   e  <  0 

The  derivation  of  solutions  for  the  negative  £  is  almost  identi- 
cal to  3.2.1.   The  equation  (3.7)  yields 

(-j^O   =  P  (a-w)  (b-w)  (w-c) 

where 

o    K|C                v  -C  Z 
2      '  o     2  .  ,,   .       o 
P   =  -3^—  ,  P  (a+b+c)  =  — . 

Obviously  bounded  solutions  exist  only  in  [c,b]  (Fig.  3.5). 
The  periodic  solution  is: 


w(z)  =  c+(b-c)sn  (r~^ Pz,k) 


2  \I  3.  —  c 
=  a-(a-c)dn  (f— - —  Pz,k) 


or 


.  2    b-c 

k   =  

a-c 


2  2K 
u(x-vt)  =  c+Asn  ( — (k   x-wt),k)  (3.14) 


=  c+  ^rfl-dn(%k   x-0)t),k)3 


where 


k2 

v  n 

,~x  -   w,-/ j 

A  =  b-c 

(ampli  tude) 

k 2. 

x   4kK  sj 

/■$* 

(wave  number) 

co  =  vk 

X 

(frequency) 

K  =  K(k) 

and 

v2  =  C  2(1-  -Lfi  (a+b+c))  .  (3.15) 

o      J 

Obviously  the  velocity  becomes  "supersonic"  if  a+b+c  <  0,  and  vice  versa. 

We  can  obtain  the  more  explicit  (two  parameter)  solutions  for 
zero  background  waves  (u=0)  by  following  the  same  procedure  as  before. 
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Fig.  3.5.   Quadratic  nonlinearity  (&  <  0) :   periodic  wave, 
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Results   are 


is 


u(x-vt)    =   -   —  {dn    (—  (kx   x-U)t),k)-   |) 

k 


3k2  K 


The    critical   modulus   k      is   also   "   .98.      The   dispersion  relation 

c 


2         _   2,     2  .4K/,.  2,    3E      ox,     4 

a-     =   C      k        -    h( — )    (k  +  — -  -2)k 
ox  vtt/v  Kx 


Hence 


A  2  2      2  4 

if   k   -►  0      and   -r-   finite,      jj      =   C      k      -4nk  and 

,2  o      x  x 

k 


if  k  -  1 


2    n  2,.^  'EJAN,   2 

jo   =  C   (1+  — )k 

o       3     x 


Solitary  Waves 

We  also  have  the  solitary  wave  solution  in  the  limit  a  ->b  (Fig.  3.6) 


u(x-vt)  =  b-Asech  (ax-|3t) 


(3.16) 


where 


A  =  b-c 


a  =  yiil 


and 


P  =  <  i  •  v 

2     2 
v   =  C   (1- 
o 


_llA 


(2b+c)) 


(3.17) 


For  the  zero-background  wave,  (b=0) 


u(x-vt)  =  -Asech  (ax-;  it  I 


(3.18) 


where 


A  =  -c  >  0 


/HE 
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pc«» 


Fig.  3.6.   Quadratic  nonlinearity  (e  <  0) :   solitary  wave 


and 

o       3 
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v2  "  C ,2(1+-TA)-  (3-19) 


Only  the  compressive  solitary  wave  is  possible  if  £  <  0.   Com- 
bining this  result  with  the  one  obtained  for  £  >  0,  we  can  summarize  that 
the  solitary  wave  is  always  created  in  the  region  in  which  the  nonlinearity 
becomes  hard .   This  is  one  of  the  common  features  of  the  solitary  waves. 

3.3   Cubic  Nonlinearity  (u-  ^  0) 

We  can  expect  that  the  higher  order  nonlinearity  will  produce 
more  variety  in  the  solutions.   We  often  put  £=0  in  this  section  for  the 
sake  of  simplicity,  because  the  quadratic  nonlinearity  does  not  change  the 
nature  of  the  solution.   That  is,  the  cubic  term  in  (3.7)  can  be  always 

absorbed  in  the  quartic  term  by  a  simple  transformation  of  variable: 

2         2       2    2 

^o   4    'Co   3    V  "Co    2 

P(w)  =  -  — w   - w   +  w   +  Dw  +  E 

x  '  6h        3h         H 

2  2 

P.C  4    C    2        2        2 

o 

+  D*  (w+  ^— )  +  E* 

l\i 

where  D1,  E'  are  new  constants. 

It  is  clear  that  the  existence  of  the  cubic  term  results  in  the 
shifting  of  the  solution  and  of  the  velocity. 

3.3.1   u-  >  0  (Hard  Nonlinearity)  ? 

2  2    ^Co 

Let  P(w)  =  -P  (w-a) (w-b) (w-c)  (w-d)  where  p   =      .   a,b,c,d 

are  the  four  real  roots  of  P(w)  =  0  and  they  satisfy  the  following 

equations:  « 

a  +  b  +  c  +  d  =  — ^  =  -  —  (3.20) 

3kp2      * 
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2   2 
v  -C  2 

ab+bc+ca+ad+bd+cd  = -j_  =  -  ~(~r  -  1).     (3.21) 

Kp2      ^  C  2 
o 


Periodic  Wave 


There  are  two  bounded  regions  in  w,  for  which  the  periodic  solu- 
tions exist.  They  are  [d,c]  and  [b,a]  (Fig.  3.7).  The  periodic  solutions 
are  expressed  in  terms  of  the  Jacobian  elliptic  functions  (see  Appendix  B) 


or 


and 


or 


where 


and 


a  c_d   2/y(a-c)(b-d) 

d  +  a sn  (-^ f-* L pz,k) 

/  \  a  -  c          z 
w  (z)  = 


1  +  £lA   jti2^/(a-c)(b-d) 


2       PZ'k) 


d  <  w   <  c 


,  ,    c-d    2  .2K,.       .  .  . 

d  +  a sn  ( — (k   x-uut)  ,k) 

a-c       tt   x 

"l(X"Vt)  "  ,  4.  <=-<>    V*a ^T^  '  (3-22) 

1  +  sn  ( — (k   x-uot),k) 

a-c      rr  x 


j  /a"bs    2  ^(a-c)  (b-d)     ,  . 
a  +  d(— )  sn  (^ ^ *-  pz,k) 

w   (z)  =  ~ b  <  w   <  a 

a-b   2  ,V (a-c)  (b-d)     ,  s 
1  +  —  sn  (^ ^ *-  pz,k) 


a  +  d(^)  sn2(%k   x-out),k) 
UlI(x"Vt)  -  ,    a-b    2,2K      


1  +  b^d"  sn  (~(kx  *-«*>& 


.2    (a-b)  (c-d)  ,   , .,   N 

k  =  (a-c) (b-d)  (modulus) 


k  = 


^°   /  M.(a-c)(b-d) 


x    4K  V      6h 

uu  =  vk. 


X 


K  =  K(k)    (the  elliptic  integral  of  the  first  kind) 


44a 


Fig.  3.7.   Cubic  nonlinearity  (m-  >  0):   periodic  waves   I 
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The    velocity   is    given  by    (3.20)    and    (3.21) 


v2   =   C      [1    -  ^(ab+bc+ca+ad+bd+cd)] 
o  6 

=   Cq2[1    -   ^-((a+b+c+d)2    -    (a2+b2+c2+d2))] 

2 
=  CQ2[1   -  |-  +  Jjj(a2+b2+c2+d2)]    .  (3.24a) 

Or,    eliminating   d    completely   yields 

v2   =  Cq2[1  +  |(a+b+c)   +^-{(a-b)2   +   (b-c)2  +    (c-a)2}]    .       (3.24b) 

If  £  =  0,  the  velocity  is  always  greater  than  C  .   For  a  nonzero  e,  the 

2 
condition  e  —  3 j-L  <  0  must  hold.   This  condition  guarantees  the  stability 

of  the  lattice;  the  interaction  force  is  always  attractive  when  the  relative 

displacement  is  large  enough. 

Solitary  Waves 

In  the  limit  b  ->  c ,  both  solutions  (3.22)  and  (3.23)  approach  to 

the  solitary  waves: 

,   b-d    .2,   n   N    b(a-d)   a(b-d)    .2, 

d+a — -  tanh(ox-pt)      N    / s   l    sech  (ax-pt) 

v*-vt>  - ,;:,  -r    -  ■  l'-\   b-d     ?     R, —  <3-25> 

1+  — r  tanh  (ax-Bt)    — sech  (ax-Bt) 

a-b  K      a-b   a-b  K 


and 


a+d^  tanh2 (ax-Bt)   ^^X   -  *g£±   sech2 (ax-Bt) 

i  ,  a~p  ,   ,  2 ,    _  ,     a-d   a-b    .  2  .    _  . 
1+  r~j"  tanh  (ax-Bt)    r— r  ~  777  sech  (ax-Bt) 


UII^>  ■  ,.l,         .2.   ~  ■  l-i       a-b    ,2 <3-26> 


where 


a  =  J   T^T(a-b)(b-d) 
and 

B  =  va  . 

Hence  the  solitary  wave  solution  in  general  is  an  infinite  series  of 

2 
sech  (ax-Bt)  even  in  this  approximation.   If  we  choose  the  special  value 
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b  =  c  =  0,  the  above  solutions  become  simple.   In  this  case  we  can  obtain 
the  zero  background  solutions.   Since  a+d  =  0, 


and 


with 


t        -\       1-tanh  (qx-Bt)         ,  ,„ 
u].(x-vt)  =  -a-  ^ u-*1-  =   -a-sech(2ax-2Bt) 

1+tanh  (ax-Bt) 


u   (x-vt)  =  a-sech  (2ax-2Bt) 


(3.28) 


and 


B  =  va 


2     2     e 
v  =  C  Z(l+  -  a 
o     J 


{?•'> 


(3.29) 


Periodic  Waves  around  the  Natural  Length  of  Lattice 

Since  P(w)  is  a  quartic  polynomial  with  the  negative  leading 
coefficient,  bounded  solutions  are  possible  even  if  P(w)  does  not  have 
four  real  zeroes.   As  a  simple  example,  we  will  consider  the  periodic  wave 
which  is  symmetric  around  the  natural  length  of  the  lattice  and  also  t=0. 
Such  a  wave  can  be  realized  by  putting 

P(w)  =  (a2-w2)(w2+b2)p2  .   (Fig.  3.8) 


The  equation  is 


.dw.      2,2   2 .  .  2  ,  2 . 
(— )   =  P  (a  -w  )(w  +b  )  , 
dz 


solution  of  which  is 


/~~2   2  2 

(z)  =  acn(V(a  +b  )  pz,k),  with  k   = 


2  u2 
a  +b 


(See  Appendix  B  )   Or,  in  terms  of  x  and  t 


,2K 


u(x-vt)  =  acn(— (k   x-^)t),k) 


(3.30) 


(3.31) 


47 


PM 


—  M> 


Fig.  3.8.   Cubic  nonlinearity  (u.  >  0) 
wave  II  (symmetric). 


periodic 
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where 


rr 

kx  =  2K 


hlC 


6h 


(a2+b2)  =  f-  C   /£-  ^ 
y    2K  o  J   6h  k 


uu  =  vk 


K  =  K(k)  . 


The  velocity  is  given  by  the  quadratic  term  in  (3.30) 


2  r   2 
v  -C 
o 


2  2 

o   .  L      I  o 

(a  -b  )  = 


6>i 


(a  -  —  +a  )  =  _,  (2-  -) 

k  k 


or 


v2  =  C  2(1+^-  (2  -  ij))  . 
o       6        ,  I 
k 

Hence  the  wave  is  "supersonic"  if  k  >  —  and  "subsonic"  if  k  <  —  .   The 


>fl 


dispersion  relation  is 


2   „2.  2,.,|  2,     ,2KN   ,  4 
a.  =C   k   (1+  t1  a  )  -  * ( — )   k 
o   x     3  *i  x 


(3.32) 


2     2  2      4 
As  a  -0,  (3.32)  reduces  to  uu   =  C  k   -  *k      ,  the  dispersion  relation  for 
'  o   x     x  r 

the  linearized  equation  of  (3.5) 

ii  =  C  2u"  +  -hi""  . 
o 

Algebraic  Solitary  Waves^l  ' 

All  the  solitary  wave  solutions  discussed  so  far  (both  quadratic 

and  cubic  case)  decay  exponentially  as  z  ►  +°°.   But  there  is  also  an  extra 

-2 
solution  for  the  cubic  nonlinearity  which  behaves  like  z   as  Z  *  +00.   Such 

a  solution  exists  if  P(w)  has  a  triple  zero  (Fig.  3.9).   Namely, 


(~)  =  p2(w-a)3(b-w)  -;  P(w)    (a  <  b) 
dz 


(3.33) 


has  a  solution 


w(z)  =  a  + 


(3.34) 


l  +  (2bpz)' 
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Fig.  3.9.   Cubic  nouLinearity  (u.  >  0)  :   algebraic  solitary  wave, 
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To  simplify  the  formulation,  we  assume  that  £=0.   Since 

P(w)  =  -  p  (w  -(3a+b)w  +3a(a+b)w  -a  (a+3b)w+a  b)  , 

we  obtain  the  following  equations  from  the  cubic  and  quadratic  terms 

2 
o-(3a+b)  =  -  -~-  =  0 

2    2 

2          ^   -Co  > 
3ap  (a+b)  =  -  f—   . 

From  the  first  equation 

b 
a-  -  j  • 

From  the  second  equation 

2 
2    r  2   ^fo_  ,2    ,2. 
V   "  Co   =  ~  (3  b  }  ' 

The  solution  (3.34)  may  be  rewritten  in  terms  of  b  alone: 

/  n      b      b 
w(z)  =  -  -  +         ~      . 

l+(2bpz) 
If  we  define  the  amplitude  of  the  wave  as  A  =  a-b 

=  !>■ 

w(z)  =  ff-1  + 


3  C  2 

,     o   ,2  2 

1-K-r A  z 

8*< 


The  velocity-amplitude  relation  is 


v2  =  C2(1  +T^A2) 


16 

The  algebraic  solitary  wave  differs  greatly  from  the  ordinary 

2 
(sech   type)  solitary  waves.   The  differences  are  summarized  as  follows: 

(1)  decay  is  much  slower  (Lorenzian) ; 

(2)  the  amplitude  is  large  and  the  condition  to  allow  such  wave   is  very 
strict  (triple  zero  required  for  P(w)). 
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(3)   the  wave  must  cover  both  the  positive  and  the  negative  w  regions. 

In  order  to  obtain  the  zero  background  wave,  we  must  introduce  the  nonzero 
e,    in  which  case  the  solution  yields 

w(z)  =  -  — 


1  +  — z z 

3  ho. 

Hence  the  amplitude  is  determined  uniquely. 

(4)   The  velocity  is  slow  compared  with  the  ordinary  solitary  wave  of  the 
same  height. 

Since  the  algebraic  solitary  wave  exists  under  very  limited 
conditions,  we  will  exclude  this  type  of  solution  in  the  future 
consideration. 

3.3.2   u.  <  0  (Soft  Nonlinearity) 

We  can  put 

2 
P(w)  =  p  (w-a)  (w-b)  (w-c)  (w-d) 

where 

2    MCo2 

a+b+c+d  =  -At  (3.36a) 

IJH  22 

2  v  -C 

p  (ab+bc+cd+da+ac+bd)  =  —  .  (3.36b) 

Periodic  Wave 

There  is  only  one  region  of  w  in  which  bounded  solutions  exist 

(Fig.  3.10).   The  solutions  are: 

a  /hz£.\      2  /V(a-c)  (b-d)     .  . 
c-d.(^)sn  Q"    '±  pz,k) 

w].(z)  =  — 


^^(a-cHb-d)  pzk) 
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PM 


Fig.  3.10.   Cubic  noniinearity  (m-  <  0):   periodic  wave, 
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or 

b-a  < 

a-c 


l   ,P-cN   2  /v/Ca-c)  (b-d)      ,  . 
b-a( )sn  (^ ^ <—   pz,k) 


wn(z)  = 


/b-c.   2  ,y(a-c)(b-d) 
1  -  (~ )sn  <~ ^ *-   Pz,k) 


As  is  discussed  in  Appendix  C,  these  two  solutions  differ  from  each  other 
only  by  phase.   That  is,  we  can  match  w  and  w   by  shifting  one  of  the 


lem 


by    2K/ps/(a-c)(b-d)     . 


Going   back   to    the   x,t   coordinates,    w      and   w        are,    respectively, 


c-d(^7)sn2(—  (k      x-uut),k) 
b-d  rrv  x 


UI(X_Vt)    =  : ,b-c,      2,2K„ ~  (3'37a) 


and 


1         (— )sn    (—  (kx   x-u)t),k) 


b-a( )sn    ( — (k     x-0)t),k) 

un(x-vt)  - ;    'b-c,  2,2k" —  <3-37b> 

1    -    ( )sn    (—  (k     x-o)t),k) 

a-c  rr    x  ' '    ' 


where 


C     TT 


_  _o_    /   u,(a-c)(b-d) 
x  "      4K  V  6h 

0)  =   vkx 

K  =  k(k)    . 

The    velocity   is    given   by    (3.36a)    and    (3.36b): 

2   „   2r.     e2     LL  2  L2   2  ,2^ 

v  =  co [1  +  iR  "  "ri(a  +b  +c  +d )] 


:Q2[1  +  j(a+b+c)  -  -^-{(a-b)2+(b-c)2+(c-a)2}].   (3.38) 


2     e2 
We  can  see  that  the  velocity  is  always  less  than  C   (1+  r-i — r)  .   In  par- 

o     3  |  p.  I 

ticular  the  velocity  is  less  than  C   if  e  =  0.   We  should  note  the  periodic 

o  r 

wave  solution  is  always  unstable  for  a  large  enough  amplitude  when  u-  is 
negative.  Therefore  the  only  amplitude  that  yields  the  right  hand  side 
of  (3.38)  positive  is  admissible. 
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Periodic  Wave  around  the  Natural  Length 

Since  we  consider  a  symmetric  wave  around  w=0,  we  put  b=-c  >  0, 

hence  a=-d  >  0  (Fig.  3.11).   We  also  assume  that  e=0  for  the  sake  of 

simplicity.   Therefore 

.  2   4ab 

k  =  2 

(a+b) 

C  TT    i 

x    4K  J   6k 


and 

=   C   2[1    -  f  {(afb)2+2bZ-2ab  }]  =   C   Z[l    -   £{(1-  \)  (aZ+bZ)+2bZ  }  ]   <  C 


2  ,r,        H-  r,  k2w  „,  2        ^  „  : 

v 

O      "  t>   "  '  O  DZ  O 


From    the   above    three    equations,    we   have 


2  -  C  2Cl  -  *(1-  t-><«=)     S5.k2.ttb2] 


o      6  v    2  VC  tr    u  x    3 

o 

.  C  2(1  -  f  b2)  -  «A  (I"  ^T)k  2,  (3.39) 

O         J  ' '  Z    X 

or  2 

2    „  2,  2    ^Co      ,4K  2  .   k2   4  „  ,  m 

ud   =  C   k    -  — —  -  h(— )  (1-  — )k    .  (3.40) 

ox      3         n       I      x 

Solitary  Wave 

If  either  a  approaches  to  b  in  (3.37a)  or  d  approaches  to  c  in 
(3.37b),  we  have  the  solitary  wave  solution  on  the  positive  (rarefactive) 

background  or  on  the  negative  (compressive)  background,  respectively. 

2       2 
Since  sn   *  tanh  ,  these  solutions  are,  respectively, 

c-de^Hanh  (ax-pt) 

u(x-vt)  =  ^ 5 (3.41a) 

1  -  (^)tanh  (0K-pt) 


and 


b-a(— )tanh2(ax-(  I  I 


u(x-vt)  =  ^ ~ (3.41b) 

1  -  (— )tanh  (ox-^t) 
a-c 


55 


Fig.  3.11.   Cubic  nonlinearity  ((J,  <  0) :   periodic 
wave  (symmetric). 
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where 


rCQ2(b-c)(b-d) 
a  =  -* ^ *-   P  =  / — for  (3.41a) 


(a-c)  (b-c)(j,C 

2^ —  for  (3.41b) 


P  =  VO! 


and 


0  0  o 

v   =  C   (L-  £(b-c)  )   for  both  cases 
o      6 

2    u  2 

Note  that  there  is  no  solitary  wave  solution  on  zero  background:   if 

a=b  were  0,  it  would  give  a+b+c+d  -    c+d=0.   But  this  is  contradictory 

to  d  <  c  *   0.   On  the  other  hand,  if  c=d  were  0,  it  would  give  a+b=0 

which  contradicts  the   condition  a  >  b  ->  0 .   Hence  the  above  proposition 

is  proved.   This  corresponds  to  the  fact  that  — ,  —  , where  f(x)  is  the 

dx 

force  function,  takes  the  maximum 
Extra  Solution 

There  is  also  an  extra  solution  for  K  ■  0.   But  this  solution 
is  quite  different  from  the  algebraic  solitary  wave  for  ,-  >  0:   it  requires 
the  positive  background  on  one  side  of  the  lattice  and  the  negative  back- 
ground on  the  other.   Since  we  are  not  interested  in  this  solution,  we 
will  give  only  a  brief  description.   We  assume  that  P(w)  has  two  double 
roots  a  and  -a  (Fig.  3.12).   The  equation  is 

2 
,dw.      2  .    .  2  .    .2 
(— )   =  P  (w+a)  (a-w) 
dz 

2   2   2  2 
=  P  (a  -w  ) 
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Fig.  3.12.   Cubic  nonlinearity  (|i  <  0) 
wave  solution. 


non-solitary 
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the  solution  which  is 

w  =  a • tanh  (apz) 


v2  =    C  2(1  -  ^  ) 
o       3« 


Table  3.1  summarizes  all  the  solutions  discussed  in  »3 . 2  and 
<3.3. 

3.4   Piecewise  Linear  Approximation  and  Explanation  of  the  Existence  of 
Solitary  Waves 

So  far  we  have  treated  the  polynomial  type  force  functions  exten- 
sively.  We  have  made  use  of  the  special  properties  of  the  elliptic  func- 
tions in  deriving  the  solutions.   Some  of  the  common  characteristics  of 
the  solutions  are: 

(1)  the  periodic  solution  can  be  either  supersonic  or  subsonic  depending 
on  the  wave  form; 

(2)  the  solitary  wave  solution  is  a  special  case  of  the  periodic  solution; 

(3)  the  solitary  wave  is  always  formed  in  such  a  way  that  the  interaction 
force  becomes  stiffer  as  the  displacement  increases; 

(4)  the  propagation  velocity  of  the  solitary  wave  is  always  higher  than 
the  ambient  sound  velocity; 

(5)  the  wave  form  of  the  solitary  wave  may  be  approximated  by  the  exponen- 
tial function  in  the  region  where  the  displacement  is  small. 

On  the  other  hand  we  learned  in  Chapter  2  that  the  harmonic  lattice  has  the 
exponential  solutions  as  well  as  the  sinusoidal  solutions.   This  fact 
together  with  (3)  and  (5)  suggests  the  lorce  described  by  two  linear  func- 
i  inns  as  a  simple  model: 


59 


x> 

H 


C 
c« 

CM 


c 

•H 
CO 

C 

o 
•l-l 

4-1 

a 
r-i 

o 

CO 


CO 
CO 

o 
a. 

<u 

4-1 


U-l 

o 

cfl 

B 
E 
3 

CO 


O 
co 

ctj 
M 

4-1 
CD 


<D 
> 

>i 

V-l 
CO 


O 
•i-l 

u 


O 


o 
II 

ri 


c73<i;Qpi<HHO 


o 


O 
II 


U    S    CQ    H    U 


60 

f  (x)  =  x         x^b 

=  (1+A)x  -  Ab   x  <  b  .    (Fig.  3.13) 

The  equation  which  corresponds  to  (3.5)  in  each  region  of  u  is  the  follow- 
ing: 

u  =  C  2u  +  Hu""  u  >  b 

o 

u  =  C  2(l+A)uM  +  m""    u  <  b  . 
o 

We  assume  that  the  solution  u(x-vt),  its  first  derivative  and  its  second 

derivative  are  continuous  at  u=b .   Putting  w(z)  =  b-u(x-vt)  and  integrating 

the  equations  for  w(z)  twice  with  respect  to  z  yield 


2    2 
hw"  =  (v  -C   )w  +  D.z  +  E.       (w   0) 
o       1     1 

HW"  =  [v  -C   (1+A)]W  +  D  z  +  E   (w  >  0) 


(3.42a) 
(3.42b) 


where  D1 ,  E  ,  D„ ,  and  E?  are  the  integration  constants.   As  we  discussed  in 
.^3.1,  we  can  put  D   =  D   =  0.   From  the  condition  that  w"  is  continuous  at 
w=0,  we  have  E.  =  E   =  E.   Multiplying  the  both  hand  sides  of  (3.42a)  and 

(3.42b)  by  w'  and  integrating  with  respect  to  z  yield 

2 
^w*2  =  (v2-C  2)t  +  Ew  + 


2* 


o  '2 


w  ■  0 


M    2     r  2     9        1W 

rw1   =  [v  -C  "(1+4)]-  +  Ew  +  F0   w  ■  0 
I  o         I  I 


(3.43a) 
(3.43b) 


where  F1  and  F_  are  the  integration  constants.   From  the  condition  that  w' 
is  continuous  at  w=0,  F   =  F_  =  F.   Dividing  the  both  hand  sides  of  the  above 


equations  by  ~  yields 


2-r  2 

2    V   o    2    9E    2F 
w'   =  w   +  — w  +  —  =  P,  (w)     w  <  0  (region  I) 

1  (3.44a) 

v2-Co2d^)   2   2E    2p  .n^    .    Tn 

w'   =  w  +  — w  +  —  -  P?(w)   w  >  0  (region  II) 

(3.44b) 


Co2A 

Obviously,  P0(w)  =  P,(w)  -  -^ — w2 
i.  I        I  •  < 


61 


1 

X) 

1 

1 
1 

JLj    I 

1 
1 

/                         ^r 

!      / 
-N    / 

I            1             > 
1       '     / 

0 

1  if 

/l 

/      1 

/         1 

/            1 

/               1 

/                  1 

Fig.  3.13.   A  typical  piecewise  linear  force 
function  (one  break  point). 
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The  solutions  of  the  differential  equation  of  the  forms  given  in 
(3.44)  are  known  to  be  either  the  sinusoidal  functions  or  the  hvperbolic 
functions,  the  exponential  function  being  included  in  the  latter  as  a 
special  case.   More  explicitly,  if  the  quadratic  polynomial  P  (w)  or  P„(w) 
has  the  positive  leading  coefficient,  the  solution  is  +cosh,  sinh  or  exp 
depending  on  whether  the  discriminant  of  the  polynomial  is  positive,  nega- 
tive or  zero,  respectively.   On  the  other  hand,  if  the  leading  coefficient 
is  negative,  the  solution  is  the  sinusoidal  function,  provided  that  the  dis- 
criminant is  positive;  otherwise  there  is  no  real  solution.   In  solving 
(3.44)  we  can  select  one  of  the  above  mentioned  functions  in  each  region, 
depending  on  the  leading  coefficient  and  the  discriminant  of  P  (w)  or 
P-(w),  then  connect  them  in  such  a  way  that  the  derivatives  of  the  selected 
functions  match  up  to  the  second  order  at  w=0 .   There  are  two  different 

cases  of  our  interest. 

2     2 
3.4.1   v    C    (Subsonic  Waves) 
o 

Both  P  (w)  and  P  (w)  have  the  negative  leading  coefficients. 
Obviously,  we  have  the  periodic  solution  in  (w  ,w  J  where  w   is  the  smaller 
zero  of  P„.  (w)  and  w   is  the  larger  zero  of  P.  (w)  .   (Fig.  3.14.)  The  solution 
may  be  obtained  bv  connecting  two  sinusoidal  waves  at  w=0 .   The  resulting 
wave  is  not  symmetric  with  respect  to  the  z-axis  due  to  the  difference 
between  P  (w)  and  P  (w)  .   In  order  to  maintain  the  wave  form  as  the  composite 
wave  propagates,  the  phase  velocity  of  the  wave  in  each  region  must  be  the 
same.   This  is  best  explained  graphically  by  using  the  dispersion  curves 
(Fig.  3.15).   Since  the  force  function  is  steeper  in  region  II,  the  disper- 
sion curve  for  this  region  is  always  located  above  that  for  region  I.   For 

a  given  phase  velocity  v  ,  we  draw  a  straight  line  in  the  x-k   plane,  which 

P  x 


63 


Fig.  3.14.   Subsonic  wave  in  piecewise  linear  lattice. 


u)=V-p&x 


Fig.  3.15.   Dispersion  relation  for  subsonic  waves 
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crosses  two  dispersion  curves  at  k   and  k  .   This  is  possible  for  a 
suitably  chosen  v  .   Then  we  put  the  solution  in  region  I  as 


w  (z)  =  A  cos  (k  z)  +  B 


(3.45a) 


and  the  solution  in  region  II  as 


wi;L(z)  =  +A2cos(k2z+0)  +  B2  . 


(3.45b) 


We  have  the  following  six  equations  to  determine  A  ,  B  ,  P,  A9,  B  ,  z 


from  (where  w  (z*)  =  w   (z*)  =  0): 

A   +  B   =  w   >  0 


A  cos(k  z  )  +  B   =  0 
B2  -  A2  =  w2    0 
A2cos(k2z'+^)  +  B2  =  0 

A  k  sin(k  z  )  =  A  k  sin(k  z  +d) 

2        *        2 
A  k   cos(k  z  )  =  A  k   cos  (k  z  +p) 


(3.46a) 
(3.46b) 
(3.46c) 
(3.46d) 
(3.46e) 
(3.46f) 


From  these  equations  we  can  eliminate  all  but  A  and  A_ ,  result- 
ing in  the  following  two  equations,  each  representing  a  hyperbola: 


l- A  : 

2    2    2    1 
w2  (k2  -kL  ) 


4  2    2   2 

K  ■;■  w2(k   -k   )n 

± :  a  i  ± 1 — 

2   2    2    2  •  ?        2 

k,  V  (L  -k,  )  ^       k/ 


1   2  v  2    1 


=  1 
(3.47a) 


4  2    2   2 

k  wi^k2  ~kl  ^ 

2   2~  2  .  2  !  Al+ 


k2  W]_  (k2  -ki 


k 


J- 


Wl  ^  2    1  ' 


1. 


(3.47b) 


We  can  solve  (3.47)  for  A.  and  A„  graphically  (Fig.  3.16). 

Since  there  is  only  one  intersection  of  the  two  graphs,  A.  and 
A_  are  determined  uniquely  for  the  given  set  of  k  ,  k. ,  w  and  w  .   We  can 
obtain  all  other  unknowns  in  a  straightforward  manner  from  (3.46).   By 
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Fig.  3.16.   Graphical  solution  of  (3.47) 
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changing  v  ,  it  is  possible  to  obtain  the  relations  among  such  quantities 

as  the  amplitude  (w  +w  )  ,  the  wave  number  and  frequency,  but  we  will  not 

go  into  the  detailed  discussion,  since  we  are  interested  in  showing  the 

existence  of  the  solutions  similar  to  those  in  §3.2  and  the  feasibility  of 

the  piecewise  linear  model  rather  than  the  numerical  calculation. 

3.4.2  C    <  v  <C   (1+A)   (Supersonic  Wave) 

In  this  case  P  (w)  has  the  positive  leading  coefficient,  whereas 

P  (w)  has  the  negative  leading  coefficient.   Both  the  periodic  wave  and  the 

solitary  wave  solutions  exist  depending  on  the  discriminant  of  P  (w) . 

Periodic  Wave 

The  periodic  wave  is  obtained  by  choosing  the  cosh  function  in 

the  region  I  and  the  cos  function  in  the  region  II,  then  connecting  them 

together  (Fig.  3.17).   The  dispersion  relations  are 

2  2      2  4 

uu     -    C     k        +   "k  (region    I) 

o      x  x 

0?    =    C  2(l+£)k   ''    -    "k   4  (region    II)     . 

O  XX 

They  are  also  shown  in  Fig.  3.18. 
Solitary  Wave 

As  the  discriminant  of  P,(w)  vanishes,  the  periodic  wave  reduces 
to  the  solitary  wave  (Fig.  3.19).   In  such  limit  the  wave  consists  of  two 
exponentially  decaying  curves  and  a  part  of  the  sinusoidal  curve,  since 
cosh(z)  reduces  to  exp(z)  or  exp(-z).   The  dispersion  curves  are  given  in 
Fig.  3.18. 

It  is  obvious  from  the  foregoing  discussion  why  the  crest  of  the 
solitary  wave  must  be  formed  in  the  region  II,  where  the  force  function  is 
steeper:   only  in  such  case  can  the  sinusoidal  part  of  the  wave  propagate 
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— z 


Fig.  3.17.   Supersonic  wave  (periodic)  in  piecewise  linear  lattice. 
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Fig.  3.18.   Dispersion  relation  for  supersonic  wave 


Fig.  3.19o   Solitary  wave  in  piecewise  linear  lattice. 
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as  fast  as  the  exponential  part.   Also  note  that  the  combination  of  the 
exponential  function  and  the  sinusoidal  function  is  possible  only  in  the 
supersonic  case. 
3.4.3   Summary 

We  have  shown  that  all  the  features  of  the  waves  in  the  anharmonic 
lattice  (in  continuum  approximation)  listed  at  the  beginning  of  this  section 
can  be  elementary  explained  by  the  piecewise  linear  force  model.  We  believe 
that  this  approach  is  quite  feasible. 

As  is  evident  from  Fig.  3.14,  Fig.  3.17  and  Fig.  3.19,  approximat- 
ing a  quadratic  force  function  by  two  linear  functions  is  equivalent  to 
approximating  a  cubic  polynomial  by  two  parabolas  in  the  P(w)-w  plane.   We 
can  obtain  a  better  approximation  to  a  given  force  function  by  extending 
this  approach  to  include  more  break  points.   For  example,  a  cubic  force 
function  in  §3.3  is  essentially  approximated  by  three  linear  functions 
(Fig.  3.20).   It  is  clear  that  the  compressive  solitary  wave  and  the  rare- 
factive  solitary  wave  exist  as  well  as  the  periodic  wave  for  such  function. 
The  only  exception  is  the  algebraic  solitary  wave,  since  it  stems  from  the 
multiple  zero  of  order  three  or  more  of  the  polynomial  P(w)  in  the  equation 

The  analytic  approach  in  §3.2  and  §3 . 3  will  be  very  difficult  for 
the  higher  order  polynomial  P(w),  since  we  are  quite  dependent  on  the  special 
properties  of  the  Jacobian  elliptic  functions.   One  of  the  advantages  of 
our  approach  in  this  section  is  that  we  can  obtain  the  approximate  solutions 
for  the  higher  order  polynomials  as  well,  since  we  are  dealing  with  elemen- 
tary functions. 

We  will  give  some  of  the  numerical  examples  in  $5.5.   We  will  con- 
sider I  he  hardsphere  and  harmonic  force  in  Chapter  5,  which  is  an  extreme 
case  of  our  approximation  here. 
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Fig.  3.20.   Solitary  waves  in  piecewise  linear  lattice  (two  break  points) 
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3.5   Toda's  Solution  for  the  Exponential  Lattice 

The  solutions  obtained  in  the  continuum  approximation  are  valid 
only  when  the  wave  form  varies  gradually.  The  only  exact  solutions  known 
for  the  anharmonic  lattice  are  those  obtained  by  Toda.  In  this 

section  we  briefly  summarize  the  results  obtained  by  Toda  and  compare  them 
with  the  solutions  obtained  in  ^3.2. 

Toda  considered  the  interaction  force  which  is  defined  as 

f(x)  =  ^(l-exp(-bx))  (3.48) 

where  b  is  a  parameter  which  determines  the  degree  of  nonlinearity .   If 
we  consider  the  lattice  with  an  infinite  number  of  atom,  the  first  term  in 
(3.48)  vanishes;  Toda's  force  function  is  essentially  the  repulsive  force. 

f(x)  may  be  expanded  in  the  Taylor  series 

2 
f(x)  =  x  -  |  x2  +  ^-  x3  -  ■••  .  (3.49) 

Hence,  exponential  force  may  be  approximated  by  the  quadratic  nonlinearity 
if  x  is  sufficiently  small. 

The  exact  equation  of  motion  is 

2 
d  u    .    -bu      -bu  , .     -bu   . 

n   1,„     n       n+1        n"l\  /q  cn\ 

=  ~(2e     -  e       -  e      )  ,  (3.50) 

dt 

where  u   =  v    -y   is  the  difference  of  the  displacement  of  two  neighboring 
n    ■  n+1   n 

atoms . 

Toda  found  two  particular  solutions  of  different  types,  one  being 
the  solitary  wave  and  the  other  being  the  periodic  wave.  The  solitary  wave 
solution  is  given  as 

u   =  -  ^log(l  +  j  2sech2(>,n  +  it))  (3.51) 

whi  i  -   i  =  sinha.  (3 .  52) 
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Note  that  the  relation  (3.52)  is  identical  to  the  one  given  for 
the  unbounded  solution  of  the  harmonic  lattices  (§2.1).   This  is  not  a 
coincidence.   Toda's  solution  (3.51)  approaches  to  the  exponential  func- 
tion as  an  +  (3t  -*   co.   That  is,  the  outskirts  of  the  solitary  wave  behave 
like  the  wave  in  the  harmonic  lattices,  for  which  B  =  sinhct  holds  exactly. 
Since  Toda's  solution  is  a  traveling  wave,  the  same  relation  between  a  and 
6  should  hold  all  over  the  wave.   Hence,  (3.52)  seems  to  hold  quite 
generally  for  the  solitary  waves  in  the  anharmonic  lattice. 

If  we  expand  (3.51)  and  (3.52)  in  the  Taylor  series,  the  first 

order  terms  are 

u  =  -R2sech2(an  +  St)      (R2  «  1) 
n 

and 

3  2  2 

S2  =  (crt-^)  ~  a2(l+^")  . 

From  the  latter  equation  the  velocity  yields 

2        2 

v     2   i.  i-     3    , 

a 

2 
which  matches  (3.13)  if  we  put  c   =1  and  £  =  1. 

As  for  the  periodic  wave,  the  exact  solution  is 

un  =  -log[l  +  (2S2)  {dn2(2^(kxniu)t)  -  |3]  (3.53) 

where  UJ  and  k  are  related  in  the  following  way: 


x 


s  n  ( —  k  ) 


2KC0 5L^  tt  ^x 

TT 


(3.54) 


7l+  (|>  l)sn2(^kx) 


The  relation  between  the  phase  velocity  and  the  wave  number  is 
obtained  from  (3.54)  and  is  shown  in  Fig.  3.21.   The  curves  are  very  close 


^ 
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Fig.  3.21.   Phase  velocity  vs  wave  number  relaf.ion  for  Toda '  s 
lattice . 
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to  that  for  the  harmonic  lattice  if  k  is  small.   As  k  increases  the 
velocity  starts  to  take  values  greater  than  C  for  the  lower  wave  numbers. 

Clearly  the  periodic  wave  for  Toda ' s  exponential  lattice  can  be  either 

? 
2  2K(J0 
supersonic  or  subsonic.   The  Taylor  expansion  is  valid  when  k  ( )   «  1 

2   E      2 
since  the  amplitude  of  (dn  -  — )  is  k  .   Therefore  either  k  or  ())  should  be 

K 

small  compared  with  1.   Then  the  first  order  approximation  is 

un  -  -(^)  [dn2{^(0)t+kxn)}  -  |]  .  (3.55) 


The   dispersion   relation   yields 


2         (2K  2    _    1*£(2K  < 

2K0J  l   _.    •    tt  V  3     '   TT  V 


,  2K  ,     ,2        .likl.E        1W_V 
W    (~kx}      "    (~l~  +  i-    1)(— > 

4 
=    (^  kx)2   +  }(2    -   k2    -  ^fX—^)  (3-56) 

where  we  have  made  use  of  the  following  approximation: 

2 

sn(x,k)  ™   x  -  — —  x         (k.x  «  1)  . 

o 

We  can  easily  show  that  (3.56)  matches  (3.11)  if  we  put  C   =  1,  h  =  — . 

It  is  expected  that  the  Taylor  expansion  of  (3.53)  up  to  the 
second  order  matches  the  solution  for  the  cubic  nonlinearity  (§3.3),  pro- 
vided that  the  latter  solution  (which  is  in  the  form  of  the  rational  func- 
tion) is  also  expanded  in  the  Taylor  series  and  the  terms  with  the  order 
of  sech   and  higher  are  truncated.   Detailed  calculation  will  not  be  made 
here . 

Toda  also  discovered  the  two  solitary  wave  solution  for  the 
exponential  lattice.   He  could  show  that  there  is  a  phase  shift  after  the 
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two  such  waves  collide  each  other.   The  phase  relation  is  such  that  if  two 
waves  make  a  head-on  collision  the  phase  shift  is  positive  for  both  waves, 
whereas  if  one  wave  takes  over  the  other,  the  phase  shift  for  the  faster 
wave  is  positive  and  the  phase  shift  for  the  slower  one  is  negative.   The 
important  thing  is  that  these  waves  preserve  their  identity  even  after  the 
collision,  unlike  the  conventional  shock  waves  in  the  dissipative  medium. 
Recently  Wadati  and  Toda  showed  that  both  the  original  equation 

of  motion  (3.46)  and  its  continuum  approximation  (3.5)  with  ^=0  give  exactly 

(32) 
the  same  phase  shift  for  the  two  wave  collision  problem.      Therefore,  we 

can  expect  that  the  continuum  approximation  gives  quite  good  quantative 

results  as  well  as  qualitative  results. 

Zabusky  discussed  in  his  recent  report  that  a  travelling  sine 

wave  in  the  anharmonic  lattice  breaks  into  many  small  solitary  waves  of 

different  amplitudes  and  velocities  ("L  Soliton").   Since  Zabusky  chose  the 

periodic  boundary  condition,  we  can  expect  that  these  small  waves  are 

actually  the  superposition  of  the  periodic  waves  of  different  amplitudes 


ind  wave  numbers,  each  of  which  are  similar  to  Toda ' s  periodic  solutions 


(33) 
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CHAPTER  4.   NUMERICAL  METHODS  AND  RELATED  TOPICS 

As  we  introduced  in  §1.2,  our  model  consists  of  two  essential 
parts:   the  crystal  lattice  described  by  a  set  of  ordinary  differential 
equations  and  the  heat  reservoirs  described  by  random  sampling  of  the 
velocities.   We  need  only  the  former  part  to  investigate  the  transmission 
properties  of  the  lattice  (deterministic  problem). 

As  for  the  lattice  we  have  chosen  the  free  boundary  conditions 


on  both  ends,  thus  obtaining  the  following  equations 


A, 


dt 

d2y. 


2  =  f  (y2-y].) 


(left  end) 


dt 

d2y, 


dt 


7=  f(y  ,,-y  )-f(y  -y  -, )    (n=2~N-i) 

2     ^n+1   n     n   n-1 


(4.1) 


=  -^vW 


(right  end) 


where  y   is  the  displacement  of  the  n1^  atom  from  its  equilibrium  position 
and  f  is  given  as 


2     3 
f  (x)  =  x  +  ax  +  px 


(4.2) 


We  can  also  write  (4.1)  in  a  vector  form: 

d2 

—7  y  =  F(y) 


dt 


where 


r> 


i      s 


yN  ! 

-J 


f(y2-yL) 


y9    |        »       F  =:     f(y,-y9)-f(y9-y,)  | 


3  J2'      w2  JV 


•^vW 
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We  also  use  the  vector  notation  in  this  chapter  when  we  mention 
the  system  of  ordinary  differential  equations  in  general,  in  which  case  F 
may  depend  on  y,t  as  well  as  on  y.   If  F  depends  on  y  only  we  call  such 
system  of  equations  "special  equations"  according  to  Henrici.       Through- 
out this  chapter  O.D.E.  stands  for  the  ordinary  differential  equations. 

In  §4.1  and  24.2  we  discuss  the  numerical  integration  schemes. 
34.3  we  discuss  the  heat  reservoirs  more  in  detail. 

In  ?4.4  we  discuss  the  computer  system  which  we  used  and  also 
the  semi-interactive  features  of  our  programs.   In  §4.5  we  show  the  typical 
structure  of  our  programs  by  giving  an  example.   (Also  in  Appendix  D.) 

4.1   Numerical  Integration  of  the  System  of  O.D.E. 

Generally  speaking,  there  are  two  different  types  of  numerical 
integration  schemes  for  the  system  of  equations  such  as  given  in  (4.1): 

(1)  One  step  methods  (Runge-Kutta  type) 

(2)  Multistep  methods  (Predictor-Corrector  Methods). 

As  for  the  special  equations,  we  can  solve  such  equations  directly  without 
reducing  them  into  a  system  of  first  order  O.D.E.  with  2N  variables.   There- 
fore we  have  four  different  types. 

First,  we  compare  the  one  step  methods  with  the  multistep  methods. 
(1)   The  one  step  methods  are  self-starting,  that  is,  we  only  have  to  specify 
the  initial  conditions.   On  the  other  hand,  multistep  methods  require 
the  information  of  several  time  steps  before  since  in  these  methods 
we  solve  the  finite  difference  equations.   Especially  in  our  problems, 
as  will  be  clarified  in  §4.3,  we  have  to  restart  the  integration 
everytime  the  end  atoms  interact  with  the  heat  reservoirs.   Therefore 
the  mu lii step  methods  are  very  difficult  to  use. 
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(2)  One  step  methods  usually  require  more  number  of  function  evaluations  . 
per  time  step  than  the  multistep  methods  do.   A  typical  fourth  order 
one  step  method  requires  four  function  evaluations  per  time  step, 
whereas  a  comparable  multistep  method  requires  only  two  function 
evaluations  per  time  step.   In  most  cases  this  is  the  major  obstacle 
against  one  step  methods,  since  the  function  evaluation  is  the  most 

time-consuming  part.   But  in  our  present  problems,  as  Waters  extensively 

(35) 
discussed,      the  function  which  we  deal  with  is  the  polynomial  (cubic 

at  most)  and  is  very  easily  evaluated.   Hence  the  one  step  methods  are 
quite  promising.   Gear  also  mentioned  this.      As  for  Toda ' s  poten- 
tial, however,  this  is  not  the  case;  if  we  need  to  do  a  long  time 
computation  with  Toda ' s  potential,  the  mixture  of  the  multistep 
method  (for  the  inner  atoms)  and  the  one  step  method  (for  those  atoms 
which  are  near  the  boundary)  will  yield  faster  computation. 

(3)  One  step  methods  can  be  programmed  very  simple. 

(4)  Error  bounds"   for  the  multistep  methods  are  more  explicitly  given 
than  those  for  the  one  step  methods.   For  the  one  step  methods  we  have 
to  evaluate  the  error  by,  for  example,  the  step  reduction  technique. 

From  the  above  consideration  we  will  solely  discuss  the  one  step 
methods  for  the  integration  of  the  equations  (4.1). 


* 

Remark  given  by  C.  W.  Gear, 


There  are  two  different  types  of  errors  in  the  numerical  integration 
of  O.D.E.   One  is  the  truncation  error  which  stems  from  the  approximation 
formula  for  integration.   The  other  is  the  round-off  error  which  stems 
from  the  finiteness  of  the  length  of  the  registers  in  the  computer.   Here, 
we  are  solely  concerned  with  the  former. 
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4.2   One  Step  Methods  for  the  Special  Equations 

Usual  technique  to  solve  the  second  order  O.D.E.  is  to  reduce 

them  to  a  system  of  first  order  O.D.E.  by  introducing  extra  variables  z  =y  . 

n  ^n 

There  are  quite  a  few  integration  schemes  for  the  system  of  first 
order  O.D.E.   Most  popular  ones  are  Runge-Kutta ' s  fourth  order  method  (to 
be  referred  to  as  RKC4  in  the  following  discussion)  and  its  variations  by 
Ralston  and  by  Gill.   These  methods  are  listed  in  Table  4.1.   Ralston's 
method  was  derived  as  a  result  of  minimizing  the  truncation  error  of  order 
h  ,  where  h  is  the  step  size  of  integration.       This  method  gives  the 
truncation  error  which  is  about  1/2  of  that  for  RKC4.   At  the  same  time, 
however,  the  simpleness  of  the  integration  scheme  is  lost  and  more  compu- 
tation time  results.   Gill's  method  was  derived  from  the  viewpoint  of 

/TON 

reducing  the  size  of  temporary  storage  rather  than  an  error  consideration. 

Blum  later  showed  that  the  same  reduction  of  temporary  storage  is  possible 

(34) 
for  RKC4  too.      Therefore  Gill's  method  is  quite  irrelevant  nowadays. 

(Note  that  the  simpleness  of  integration  scheme  is  also  lost  in  Gill's 

method.)    Hence  we  can  naturally  select  RKC4  as  the  most  optimal  one  among 

* 

the  fourth  order  one  step  methods.     For  a  higher  precision  computation 

Waters     and  Pasta  et  al_.      recommend  Butcher's  fifth  order  method.   It 

(18 ) 
is  reported      that  the  truncation  error  could  be  reduced  as  small  as  the 

round  off  error  with  a  choice  of  h  =  .01.   This  is  a  "super  accurate" 

example.   Usually  the  truncation  error  is  several  orders  larger  than  the 

roundoff  error.   We  do  not  require  such  a  high  accuracy  in  our  computation 


According  to  (39),  p.  363,  the  Runge  Kutta  subroutines  in  most  computer 
libraries  still  employ  Gill's  method.   Probably  this  is  just  because  the 
original  (and  the  best!)  method  was  derived  in  1901,  whereas  Gill's  method 
was  derived  in  1951,  the  latter  giving  an  impression  that  it  is  quite 
i  mproved . 
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Classical   Runge-Kutta    (1901) 


Table   4.1 
(34) 


W    -   *n  +  t(«l   +  2t2   +  2lt3   +  *4> 


K.    =    F(y  +  hK„) 
4  n  3 

(37) 
Runge-Kutta-Ralston  (1965) 


y    ,      =  y  +  h(.17476028K1    -    .55148053   K_ 
•'n+l  n  1  2 


+  1.20553547K-   +   .17118478K. ) 
3  4 


KL    =   F(yn) 


K_    =   F(y  +    .4hK.  ) 
z  n  1 


K  =  F(y  +  ,29697760hK„  +  .15875966hK  ) 

3  Jn  2  2 

K.  =  F(y  +  .218100381^.  -  3 .05096470hKo  +  3 .83286432hico) 

4  n  1  2  3 


Runge-Kutta-Gill  (1951) ^34^ 


?n+l  =  V  |(^+(2W2)K2  +  (2-V2)K3  +  K4) 

*2  =^n+2KV 

K3  =  F(yn+  IC-l-K/2)^  +  |(2-V/2)K2) 


Table  4.1  (continued) 


Nystrom's  method  for  the  special  equation  (1925) 

,  2 


(34) 


Vfl    =  yn  +  h'K  +  I^(23K1   +  75K2    ~    27K3   +  "V 
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?Ul   =   K  +  I9T(23K1   +   125K2    "    8LK3   +  125K4} 

h   =   f^n} 

K2    =    F(yn+  fh?;      +  ^  h2KL) 

h  =  ?<?n+  l^n     +     I  h^l} 

*4  -  ?<yn+ ftf;   +^h2(K1+K2)) 


Butcher's  fifth  order  method  (1964) 


(35) 


■   y  +  ^(7K,    +  32K,   +   12K,    +  32K_   +  7K,) 


n+1        Jn      90 


KL    =    F(yn) 


K2    =   F(yn-   ?1) 


?3  =  ?(vn(  5KrK2)) 


K      =    F(y   +  7    (-3K.    +  K0    +  4K_)) 
4  n      4  1  2  J 


3h 


K5   =    F(^n+I?(K1    +  V} 


*6   =   f(yn+      7(~*2   +   Uh    -    Uh   +  8^5)} 


therefore  a  good  compromise  is  the  RKC4  or  its  equivalence  in  the  sense 
of  computational  complexity  (that  is,  the  number  of  function  evaluation 
per  time  step,  the  size  of  temporary  storage,  etc).   Butcher's  fifth  order 
method  is  also  given  in  Table  4.1. 

An  alternative  approach  is  to  integrate  the  second  order  O.D.E. 
without  reducing  it  to  the  first  order  (direct  method) .   Henrici  discussed 
that  the  order  of  the  truncation  error  for  y  do  not  match  that  for  y  in  the 
general  cases  in  which  F  depends  on  y  as  well  as  y.       Therefore  the 

direct  methods  are  not  recommended.   But  things  are  different  for  the 

(34) 
special  equations.   Then,  as  Henrici  pointed  out,      not  only  the  trunca- 

— ►     ~~ * 
tion  error  becomes  the  same  order  for  y  and  y  but  also  it  is  of  the  higher 

order  error  than  the  corresponding  Runge  Kutta  methods.   For  example,  two 

function  evaluations  yield  a  third  order  method,  three  function  evaluations 

yield  a  fourth  order  method  and  four  function  evaluations  yield  a  fifth 

order  method.   We  are  particularly  interested  in  the  fifth  order  method 

due  to  Nystrom,     which  is  given  in  Table  4.1.   RKC4  and  the  Nystrom's 

fifth  order  method  for  the  special  equations  (to  be  called  NYS5)  may  be 

compared  as  follows: 

(1)  NYS5  yields  higher  order  global  error.   That  is,  the  total  error  at 
time  t  due  to  the  truncation  behaves  approximately  as  C  • h   for  a  given 
integration  step  h,  where  C.  is  a  constant  which  depends  on  t  and  on  other 

factors,  and  has  to  be  evaluated  for  individual  problems.   As  for  RKC4  the 

4 
total  error  at  time  t  is  approximately  C_ ■ h  .   Hence  if  C.  turns  out  to  be 

comparable  to  C   for  a  given  problem  then  NYS5  yields  better  results. 

(2)  The  size  of  the  temporary  storage  is  smaller  for  NYS5.   If  we  assume 
the  nearest  neighbor  interactions  only,  NYS5  needs  three  temporary  storage 


Table  4.2 
Comparison  of  NYS5  and  RKC4 
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NYS5 


RKC4 


V,  =  z 


Ki  ■  F(V 


V2  "  \   +  IK1 

«2  "  ?(*n  +  ¥2> 


V3  "  Zn  +  3*1 

*3  "  ?(7n  +  ¥3' 


V4  "  Zn  +  5(K1+V 

-►    "*,->   4h-*  „ 
K4  "  '<V  ~V 


K,  = 


V.  = 


K„  = 


F(yn) 


n    2  1 


F(yn  +  -VL) 


V-  =  z   +  -K. 
3     n    2  2 


k3=F(yn+IV2) 


V.  -  z   +  hK_ 
4     n      3 


K4  =  F(yn+hV3) 


'n+1 


n+1 


z  +  7^r(23K1+125K0-81KQ 
n   iyz     L      Z     J 

+12  5K. ) 

4 


"  yn+hZn+l"  T9l(50K2-54K3 


z   ,  =  2  +  7(K1+2K_+2K„+K.  ) 
n+1     n   6   1    2    3   4 


y  i  =  y  +  z(v1+2v0+2v„+v.) 

^n+1     n    6   1    2    3   4 


Note 


+100K. ) 
4 


V„,  V...  V.  are  not  neces' 
2    3    4 


sary  for  actual  computation, 
They  are  shown  only  to  make 


comparison  with  RKC4 

are  the  temporary  vectors 


K,,  K_ ,  K^ 


Note:   V2,  K2 ,  V3 ,  K3 ,  V4 ,  K4  are 

not  necessary  for  actual  computation 
since  we  can  accumulate  the  sums 
(K  +2K2+2K3+K4)  and  (V^V^V^) 

in  K1  and  V,,  respectively.   On  the 

other  hand  we  need  temporary  vectors, 
one  for  y*  and  the  other  for  z*.  Hence 
4  temporary  vectors  are  necessary. 
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vectors,  whereas  RKC4  needs  four.   In  general  they  are  four  and  eight, 
respectively.   As  a  matter  of  fact,  this  difference  is  not  so  critical. 
(3)   The  number  of  arithmetic  operations  (multiplication,  addition)  is 
almost  the  same  for  both  methods,  RKC4  being  slightly  superior  to  NYS5. 

To  compare  the  total  error  due  to  the  truncation,  we  have  run  a 

2        3 
test  program  with  50  atoms,  f  =  x-lOx  +  66. 7x  .   The  initial  velocities 

were  given  to  two  end  atoms  only.   The  error  in  the  total  energy  at  time 
100  is  given  in  Fig.  4.1  for  h  in  [.025,  .2].   As  is  clear  from  this  figure, 
NYS5  was  always  superior  to  RKC4  in  the  given  range  of  h.   Any  effect  of 
the  roundoff  error  was  not  observed  even  when  h  =  .025.   From  this  experi- 
ment we  have  decided  to  choose  NYS5  over  RKC4  for  our  computation. 


4.3   Heat  Reservoirs 

As  we  introduced  in  §1.2,  our  heat  reservoirs  are  defined 
dynamically,  that  is,  they  consist  of  gas  molecules  in  thermal  equilibrium. 
Take  the  left  end,  for  example.   At  each  step  of  integration  (or  at  several 
steps  as  we  will  discuss  later),  the  coordinate  y1 (t)  and  the  velocity  y,  (t) 
are  examined  to  see  if  it  will  cross  the  boundary  during  the  next  time 
interval  h.   This  is  done  by  simple  linear  interpolation  y.  (t)  +  tvy..(t). 
If  this  condition  is  met,  the  velocity  is  replaced  by  a  new  velocity, 

sampled  at  random  from  the  left  heat  reservoir  distribution  of  velocity, 

2 
_v__ 

~2T 
P(v)  =  v/T*e     .   A  similar  interaction  occurs  at  the  right  end  too.   An 

interaction  occurs  only  at  the  end  of  an  integration  step,  hence  no  problem 

of  the  discontinuities;  essentially  we  restart  the  integration  with  the  new 

initial  values  at  the  end  atoms  everytime  they  (or  one  of  them)  hit  the 

boundaries.   See  Gear.      Although  we  do  not  introduce  any  viscosity 

explicitly  in  this  model,  it  essentially  exists  dynamically  at  the  boundaries, 
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Fig.  4.1.   Comparison  of  NYS5  and  RKC4  (truncation  errors) 


This  can  be  best  understood  by  considering  two  zero  temperature  heat 
reservoirs;  the  end  atoms  lose  their  energy  everytime  they  interact  with 
the  heat  reservoirs  and  the  Lattice  will  be  cooled  off  eventually.   Some 
computational  results  are  given  in  §6.1.   We  found  that  a  more  efficient 
energy  transfer  was  possible  by  letting  the  end  atoms  move  freely  for 
several  integration  steps,  that  is,  boundary  crossing  of  the  end  atoms  is 
examined  at  several  integration  steps.   (Also  refer  to  §5.3.)    It  is  also 
possible  to  obtain  a  better  contact  by  moving  the  heat  reservoirs  inward 
(applied  pressure)  so  that  more  number  of  collisions  takes  place.   JPW 
reported  that  this  increased  the  number  of  collisions.   We  should  note, 
however,  that  this  inevitably  changes  the  degree  of  nonlinearity  of  the 
lattice,  since  the  average  interatomic  spacing  is  reduced.   In  our  computa- 
tion we  did  not  apply  any  pressure  to  the  lattice. 

Nakazawa  took  a  different  approach  by  considering  the  Brownian 
motion  of  the  end  atoms  directly.   Namely  he  introduced  the  viscosity  to 
the  end  atoms  besides  the  random  force  and  the  nearest  neighbor  interaction 
Hence  the  equation  of  motion  is  changed  to: 

d  yL  dYL 

2  =  f(y2-y].)  -  ^"dt  +  B  (t)      (left  end) 

dt 

— 2   "   f <Wn>    "    f<VVl>  (n=2   ~  "-15 

at 

^N  „,  .  % 


dt 


2  =   "f(yN"yN-l}  "    ^N~dT+  BN(0        (ri§ht   end> 


where 


and 


B1(t)B1(f  )    =   2T16(t-t') 


Vt)BN(t,)  =  2V(t-t,} 
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Nagazawa  chose  |i  and  |i   that  yield  the  optimal  contact  between 
the  lattice  and  the  heat  reservoirs  for  the  highest  mode  of  the  lattice.    ' 

4.4   Computer  Facilities  and  Programs 

The  computer  which  we  used  for  our  numerical  experiments  is 
Burroughs  B6700  system  located  at  the  Center  for  Advanced  Computation  at 
the  University  of  Illinois.   The  word  length  is  48  bits  with  39  bit 
mantissa.   Its  characteristics  are  implied  by  the  following  figures:  addi- 
tion time  .20o.s,  multiplication  time  2.00jj,s  and  division  time  10.8|j.s. 
B6700  has  a  special  architecture  designed  suitably  for  processing  ALGOL 
language.   It  is  not  primarily  designed  for  high  speed  scientific  computa- 
tions, but  its  flexible  operating  system  (MCP)  allows  users  to  manipulate 
remote  terminals,  on-line  monitoring,  etc.,  without  requiring  much  exper- 
iences in  system  programming  on  users'  part.   Fig.  4.2  shows  the  B6700 
system  as  far  as  our  computation  is  concerned. 

Our  programs  are  stored  on  a  9  track  magnetic  tape  and  are  easily 
read  on  the  disc  memory.   Editing  of  programs  can  be  done  from  the  remote 
terminals,  from  which  an  immediate  access  to  the  disc  memory  is  possible. 
The  remote  terminal  which  we  used  is  the  Hazeltine  2000.   It  is  a  character 
display  terminal  (27  rows  x  68  columns),  having  the  cursor  addressing 
feature.   This  latter  feature  allows  users  to  write  programs  which  plot 
curves  (although  very  crude  since  its  resolution  is  25x60  points)  while  the 
main  program  is  being  executed.   This  is  very  convenient  when  we  monitor 
what  is  going  on  in  the  program. 

By  making  use  of  the  terminal  just  described,  we  wrote  our  pro- 
grams in  a  semi-interactive  way.   Some  of  the  features  of  our  programs  are: 
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REMOTE  TERMINAL 


LINE  PRINTER 


Fig.  4.2.   Computer  facilities 
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(1)  A  command  to  execute  a  program  is  given  from  a  terminal.   This  terminal 
will  be  tied  up  to  the  program  until  the  program  is  terminated. 

(2)  Input  parameters  are  fed  from  the  terminal.   There  are  also  some 
options  for  printing  only  the  desired  variables  such  as  the  coor- 
dinates, velocities,  temperature  and  flux. 

(3)  Intermediate  results  are  displayed  on  the  terminal  at  an  option.   Time 
interval  can  also  be  set  at  which  we  require  displays. 

(4)  Programs  can  be  terminated  at  any  time  by  sending  a  command  from  the 
terminal  if  the  user  does  not  need  it  any  more. 

(5)  An  interactive  use  of  the  terminal  is  possible  only  when  the  program 
reaches  to  the  points  at  which  it  requires  some  information  from  the 
user.   Therefore  our  programs  are  not  interactive  in  the  strict  sense. 
But  even  semi-interactive  programs  of  this  level  still  enable  us  to 
use  the  computer  system  quite  efficiently  for  the  heuristic  research. 

Instead  of  creating  a  huge,  complicated  program  that  handles 
everything,  we  have  made  separate  programs  for  considerably  different  tasks. 
Table  4.3  gives  a  list  of  programs  which  we  have  created. 

4.5   Description  of  a  Typical  Program 

As  a  typical  example  of  our  programs  we  describe  SOL/XXIV  that 
was  written  for  computing  the  heat  flux  and  the  temperature  distribution. 
A  simple  flow  chart  is  given  in  Fig.  4.3.   In  this  program  the  heat  flux 
and  the  temperature  are  computed  and  printed  (also  on  the  terminal)  at 
every  50  seconds.   The  temperature  gradient  is  computed  by  the  linear  least- 
square  interpolation.   The  entire  program  is  in  Appendix  D. 
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Table  4.3 
List  of  Programs 


SOL/XXIV 

ZABUSKY 

IMPACT 

IMPACT2 

DISORDER 

DISORDER2 

DALT 

BRKLIN 

RESPONSE 

vise 

RKC4 

NYS5 

SPHERE 

FPU 

FPU2 

HAZELTINE 

FPU3 


heat  conduction  (JPW) 

periodic  boundary  condition  (Zabusky) 

creation  of  solitary  wave  (JPW) 

creation  of  tail  (JPW) 

one  impurity  at  n=25  (harmonic) 

impurities  at  n=25~50  (harmonic) 

impurities  alternatively  (harmonic) 

piecewise  linear  force 

cooling  process  of  the  lattice  (JPW) 

viscous  ends  (Jackson's  problem) 

Runge  Kutta  Classical  test 

Nystrom  test 

hard  sphere  •+-  harmonic 

recurrence  phenomena  (cubic  force) 

recurrence  phenomena  (quadratic  force) 

test  program  for  the  terminal 

recurrence  phenomena  (free  ends) 


INITIAL 
CONDITIONS, 
PARAMETERS 
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INTERACTION 
IN  LEFT 
RESERVOIR 


INTERACTION 
IN  RIGHT 
RESERVOIR 


NYS5 
INTEGRATOR 
(ONE  STEP) 


TEMPERATURE 


TERMINAL 


TOTAL 

ENERGY , 

FLUX 


TEMPERATURE 

LEAST 

SQUARE 


YES 


PRINT 
I 


* 


Fig.  4.3 
Flow  chart  of  a  typical  program 
(SOL/XXIV) 


END 
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CHAPTER  5.   DYNAMIC  PROPERTIES  OF  ANHARMONIC  LATTICES 

This  chapter  deals  with  the  wave  motions  in  the  anharmonic  lattices. 

Specifically  JPW  lattice  (which  is  characterized  by  the  force  function 

2      3 
f(x)  =  x-lOx  +66. 7x  )  is  considered  in  full,  but  the  piecewise-linear 

force  model  and  the  hard  sphere  plus  linear  force  models  are  also  treated. 
Since  no  analytical  solutions  are  known,  the  numerical  methods  are  employed 
and  the  results  are  compared  with  the  analytical  solutions  in  the  con- 
tinuum approximation  whenever  possible. 

The  solutions  in  Chapter  3  were  considered  in  the  special  cases;  the 
solitary  waves  and  the  periodic  waves  in  the  infinitely  long  lattices,  or 

equivalently ,  under  the  periodic  boundary  condition.   Most  of  the  numerical 

(33) 
experiments  have  been  performed  along  this  line.       In  the  heat  conduc- 
tion problem  there  are  new  factors.   Firstly,  the  lattice  is  necessarily 
bounded  by  the  heat  reservoirs.   Therefore  the  effect  of  the  boundaries 
must  be  taken  into  account.   Secondly,  the  initial  condition  is  different 
from  that  for  a  pure  solitary  wave  solution;  we  assume  the  impulse  type 
interaction  between  the  lattice  and  the  heat  reservoirs.   In  fact,  an 
impulse  produces  an  oscillatory  wave  packet  together  with  a  solitary  wave. 
The  main  purpose  of  this  chapter  is  to  clarify  the  basic  energy  transport 
properties  of  the  anharmonic  lattices  from  the  numerical  experiments.   Most 
of  the  results  in  this  chapter  have  not  been  reported. 

From  §5.1  to  §5.4  we  consider  solely  the  JPW  lattice.   In   §5.1  we 
consider  the  creation  of  the  solitary  wave  by  the  impulse.   In  §5.2  we 
will  compare  various  boundary  conditions  and  the  reflection  of  the  soli- 
tary waves  therein.   We  will  show  that  the  fixed  boundary  is  a  perfect 
reflector  but  that  the  free  boundary  destroys  the  solitary  waves.   In  §5.3 
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we  will  consider  the  two  wave  collision  process.   We  discovered  that  the 
solitary  wave  slightly  loses  its  energy  in  the  collision  process.   We 
will  also  mention  the  wave  propagation  through  the  disturbances.   In  §5.4 
we  will  discuss  the  behavior  of  the  time  averaged  kinetic  energy  of  the 
atoms  when  the  waves  propagate  through  the  lattice.   This  quantity  defines 
the  local  temperature  in  the  thermal  situation.   In  §5.5  and  in  §5.6  we 
will  consider  the  piecewise  linear  force  function  and  the  hard  sphere  type 
force  function,  respectively.   We  will  numerically  show  that  the  solitary 
wave  can  exist  in  these  models. 

The  results  obtained  in  this  chapter  will  be  used  in  Chapter  6  in 
understanding  the  heat  conduction  phenomena. 

We  should  note  that  we  will  employ  the  coordinate  space  which  is 
defined  by  the  atom  number  n  and  the  time  t  rather  than  the  spatial  coor- 
dinate x  and  the  time  t.   The  former  definition  corresponds  to  the  Lagrangian 
coordinate  system  and  the  latter  to  the  Eulerian  coordinate  system  in 
hvdrodynami .s .   The  advantage  of  the  former  is  that  the  atom  number  n  is 

h  h 

fixed  all  the  time  whereas  the  n    atom  itself  moves  in  the  one  dimensional 
space.   This  difference  is  very  important  when  we  consider  the  collision 
process.   The  approach  we  took  in  deriving  the  partial  differential  equation 
in  Chapter  3  relates  the  (n,t)  space  in  the  lattice  with  the  (x,t)  space 
in  the  continuum  approximation. 

5.1   Basic  Dynamic  Properties  of  JPW  Lattice 

5.1.1   Creation  of  the  Solitary  Wave  by  the  Impulse 

We  discovered  that  an  impulse  applied  to  the  end  atom  of  the  lattice 
produces  a  solitary  wave  and  a  tail.   By  the  latter  we  mean  a  small  wave 
packet  which  consists  of  high  frequency  waves.   A  typical  profile  of  such 
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wave  is  given  in  Fig.  5.1a.   Separation  of  the  two  parts  was  observed 

at  an  early  time  for  the  initial  velocity  greater  than  .08.   If  the  initial 

velocity  is  small  (<. 05)  separation  was  not  observed  within  the  limit  of 

our  lattice  size  (N=100) .   Such  an  example  is  also  given  in  Fig.  5.1b. 

In  this  case  the  entire  wave  propagates  as  a  (dispersive)  wave  packet. 

Presumably  separation  will  take  place  after  the  wave  propagates  for  a  long 

distance.   Zabusky  also  found  a  similar  type  of  the  wave  as  a  solution  of 

ft-1) 

the  KdV  equation.      Unfortunately,  there  is  no  analytical  solution  known 

for  such  composite  (solitary  wave  plus  tail)  waves. 

First  we  look  at  the  solitary  wave  part.   We  performed  numerical  experi- 
ments with  different  initial  velocities.   Results  are  in  Table  5.1.   Fig. 
5.2  shows  the  relation  between  the  amplitude  of  the  wave  and  the  propagation 
velocity  v.   The  solid  curve  is  the  theoretical  prediction  from  the  con- 
tinuum approximation  (3.29): 

2     2     e    u,  2 
v=C(l+-A+  fk) 
o       Jo 

2 
with  C   =  1,  E  =  2x10,  M-=3X66.7  and  a=-A.   Experimental  data  are  always 

located  below  the  theoretical  curve.   This  tendency  becomes  manifest  as  the 

amplitude  increases.   Apparently  the  truncated  terms  in  the  continuum 

approximation  are  responsible  for  this  effect,  since  the  wave  form  becomes 

sharper  as  the  amplitude  increases.   Fig.  5.3  shows  the  relation  between 

the  initial  velocity  V.  and  the  propagation  velocity  v.   We  can  derive  an 

empirical  formula  v  =  1+1. 85V.  in  the  range  [.1,  .4"]. 

5.1.2   Energy  Concentration 

Since  the  propagation  velocity  of  the  solitary  wave  is  much  faster  than 

that  of  the  tail,  we  can  calculate  the  energy  contained  in  each  part 

separately.   The  energy  of  the  n   atom  is  defined  as 


Table  5.1.   Dynamic  Properties  of  Solitary  Waves 
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V. 

1 

Amp . 

V 

E   . 
sol 

tail 

Eonl  CD 
sol 

h 

Error 

.08 

.0364 

1.13 

.00292 

.00028 

91.1 

.02 

.5 

10"8 

.1 

.0465 

1.19 

.00467 

.00033 

93.4 

.05 

.9 

10"5 

.2 

.0900 

1.39 

.01960 

.00040 

98.0 

.025 

<.  5 

ID"8 

.3 

.1252 

1.57 

.04455 

.00045 

99.0 

.05 

.2 

lO"6 

.4 

.1550 

1.75 

.07952 

.00048 

99.4 

.05 

.6 

lO"6 

.5 

.1819 

1.89 

.12450 

.00050 

99.6 

.02  5 

.8 

lO"7 

.6 

.2066 

2.02 

.17940 

.00060 

99.67 

.05 

.3 

l0-5 

.8 

.2482 

2.27 

.31931 

.00069 

99.8 

.05 

.15 

lO"4 

1.0 

.2840 

2.51 

.49920 

.00080 

99.99 

.025 

.1 

.1 

.6 

lO"5 

1.  2 

.3199 

2  .  67 

1 .  5 

.3640 

2.  97 

V.:   initial  velocity 
l 

Amp:   amplitude  of  the  wave 

V:   propagation  velocity  of  the  wave 

E   , :   total  energy  in  the  solitary  wave 
sol 

E   .  •   total  energy  in  the  tail 
tail 

E,  .(%):   E   /(.5  V  2) 
sol        sol       1 

h:   integration  step 

Error:   numerical  error  in  total  energy 

means  that  such  quantity  was  not  computed. 
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Fig.  5.1.   Typical  profiles  of  waves  created  by  impulse 
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Fig.  5.2.   Velocity  vs  amplitude  relation. 
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Fig.  5.3.   Initial  velocity  vs  propagation  velocity, 
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En  =  ly/  +  *[v(yn+L-yn)  +  V^-y^)]  (5.1) 

where  V(y    -y  )  is  the  potential  energy  stored  in  the  spring  connecting 
the  n    and  n+1  st  atoms.   Results  are  in  Fig.  5.4.   Numerical  data  for  the 
JPW  lattice  are  also  given  in  the  fourth  column  of  Table  5.1.   The  stronger 
the  nonlinearity  is,  or  the  larger  the  amplitude  is,  the  more  fraction  of 
the  energy  is  concentrated  in  the  solitary  wave,  although  the  energy  in  the 
tail  itself  also  increases.   Fig.  5.5  is  a  typical  example  of  the  energy 
distribution  in  the  JPW  lattice.   Fig.  5.6  shows  the  relation  between  the 

energy  and  the  propagation  velocity.   From  this  graph  we  can  obtain  an 

2 
empirical  formula  which  is  similar  to  E  =  mv  /2 : 

E   =  M  (v-l)2/2  (5.2) 

s     s 

2 
where  M  =  .24  (L+.35(v-l)  )   (effective  mass!). 

It  is  obvious  from  the  above  results  that  the  anharmonic  lattice  can 
transmit  energy  more  efficiently  than  the  harmonic  lattice  because  of  the 
nondispersiveness  and  the  high  propagation  velocity  of  the  waves  therein. 
5.1.3   Tail  Part 

As  we  mentioned  in  5.1.1  an  oscillatory  wave  is  always  created  with  a 
solitary  wave.   Its  amplitude  is  very  small  compared  with  that  of  the  soli- 
tary wave  which  it  accompanies  with.   We  can  assume  that  the  tail  is  the 
harmonic  dispersive  wave.   This  was  verified  by  the  following  numerical 
experiment : 

(1)  Create  a  solitary  wave  and  a  tail  with  an  impulse  and  let 
them  propagate  until  they  are  separated  from  each  other; 

(2)  erase  the  solitary  wave; 

(3)  reverse  the  direction  of  the  time  and  perform  the  numerical 
integration  backward  until  t=0; 
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Fig.  5.5.   Typical  energy  distribution  in  solitary  wave, 
V.  =  .2,  t  =  37,  h  =  .05,  AE  =.5xl0-7. 


100 


Eherjy  »•»  Society  U/ave 


.01 


ool 


T 1     ■    I   I 


/ 


I      I     I    I   I  I  I 


-i — i— i — r-r 


t 

/ 


/ 


/ 


j 1 — i    i   i  i  i  | 


5oli'<J  Carve 
( Empirical  Formula.) 


Z   . 


J 1 — i  i  I  I  i  I L 


10 


V-l 


Fig.  5.6.   Energy  vs  propagation  velocity  relation  for 
solitary  waves. 
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(4)  change  the  force  function  to  the  linear  form  and  reverse  the 
time  again.   The  tail  will  propagate  toward  right; 

(5)  compare  the  wave  form  of  the  harmonic  tail  with  the  form  of  the 
anharmonic  tail  at  the  corresponding  time. 

Fig.  5.7  shows  the  result;  there  is  no  appreciable  difference  between 
them.   The  phase  delay  is  due  to  the  different  initial  condition  and  not 
due  to  the  different  propagation  velocities. 

The  tail  propagates  very  slowly  since  it  consists  of  high  frequency 
components.   It  is  quite  suggestive  that  the  interaction  between  the 
lattice  and  the  heat  reservoir  always  produces  two  waves  of  completely 
different  characters;  we  depend  on  this  fact  when  we  consider  the  origin  of 
the  temperature  gradient  (Chapter  6). 

5.1.4  Mode  Spectrum  of  Solitary  Waves 

Since  the  solitary  wave  is  formed  due  to  the  coupling  of  the  modes, 
it  is  interesting  to  see  how  the  mode  spectrum  changes  as  the  wave  propa- 
gates.  An  impulse  applied  to  the  end  of  the  lattice  gives  the  initial 
energy  distribution  as  is  shown  in  Fig.  2.2.   The  time  evolution  of  the 
spectrum  is  shown  in  Fig.  5.8.   The  lattice  consists  of  25  atoms  with  free 
boundaries.   The  initial  velocity  is  .5. 

5.2   Boundary  Effects 

We  examined  several  types  of  boundary  conditions.   The  solitary  waves 
were  injected  from  the  left  for  all  experiments. 
5.2.1   Fixed  Boundary  (Fig.  5.9) 

Solitary  waves  are  reflected  completely.   There  is  no  appreciable  energy 
left  at  the  boundary,  but  the  phase  of  the  wave  is  somewhat  advanced.   (See 
§5.3.) 
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Fig.    5.8.      Mode    spectra    (N=25)    of   solitary  waves 
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5.2.2  Free  Boundary  (Fig.  5.10) 

The  solitary  waves  decay  at  the  free  end  when  the  amplitudes  are 
small,  whereas  those  with  large  amplitudes  can  be  reflected.   This  difference 
may  be  explained  in  the  following  way.   When  the  solitary  wave  reaches  the 
end  of  the  lattice,  atoms  near  the  end  tend  to  travel  in  the  same  direction, 
resulting  in  the  expansion  of  the  lattice  there.   This  newly  created  expan- 
sion wave  can  propagate  through  the  lattice  back  again  as  a  solitary  wave 
only  when  the  force  (in  the  positive  region  of  the  displacement)  becomes 
stiff  as  the  amplitude  increases.   In  the  case  of  the  JPW  lattice,  however, 
the  force  at  first  becomes  limp  due  to  the  negative  quadratic  nonlinearity , 
then  it  becomes  stiff  again  due  to  the  cubic  nonlinearity  which  dominates 
the  quadratic  one.   Hence  there  is  a  threshold  of  the  amplitude  below  which 
the  solitary  waves  become  unstable.   To  make  this  point  clear  we  also  per- 
formed a  numerical  experiment  with  the  negative  impulses  applied  to  the 
JPW  lattice.   Fig.  5.11  shows  the  responses  of  the  lattice  to  various 
inputs.   The  critical  velocity  turned  out  to  be  V.  =  .3. 

5.2.3  Heat  Reservoir  Type  Boundary 

In  this  boundary  condition,  the  coordinate  of  the  end  atom  is  checked 
at  every  m  integration  steps.   If  the  end  atom  is  found  to  be  inside  the 
heat  reservoir  a  new  velocity  is  assigned.   We  performed  the  numerical 
experiments  with  various  m.   Since  we  are  concerned  with  the  boundary  effect 
in  this  chapter  rather  than  the  thermal  properties,  we  chose  the  predeter- 
mined values  for  V.  and  V  ,  where  V  is  the  velocity  with  which  the  end 

l      r         r 

atom  is  returned.   Specifically,  V.=l  and  V  =.1.   Results  are  in  Fig.  5.12 
a-e.   The  amplitude  of  the  reflected  wave  decreases  as  m  increases.   If  m=l 
the  end  atom  is  effectively  frozen  at  the  boundary,  therefore  the  boundary 
condition  is  very  similar  to  the  fixed  boundary  condition.   As  m  increases, 
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the  end  atom  can  receive  more  energy  from  its  left  neighbor  before  it  inter- 
acts with  the  heat  reservoir,  hence  more  energy  is  exchanged.   We  also 
examined  the  boundary  condition  (called  "soliton  eater")  which  supposedly 
yield  the  most  efficient  energy  exchange  for  the  solitary  waves;  the  end 
atom  does  not  interact  with  the  reservoir  until  its  velocity  reaches  the 
maximum  value.   The  result  is  shown  in  Fig.  5.12e.   The  effect  of  the 
boundary  condition  on  the  heat  flux  will  be  discussed  in  Chapter  6. 

5.3   Collision  Process  of  Solitary  Waves 

5.3.1   Phase  Shift  in  Two  Wave  Collision  Process 

It  was  theoretically  shown  that  the  phase  shift  occurs  in  the  two 
wave  collision  process  in  Toda's  lattice  (23.5).   In  the  JPW  lattice  the 
phase  shift  was  observed  experimentally;  both  waves  were  advanced  after  a 
head-on  collision,  whereas  the  faster  wave  was  advanced  and  the  slower 
wave  was  retarded  after  a  take-over  collision  (Fig.  5.13).   The  amount  of 
the  phase  shift  was  larger  for  the  small  amplitude  wave  in  both  types  of 
collision.   This  seems  to  indicate  that  the  motion  of  "the  center  of  the 

mass"  is  not  affected  by  either  type  of  collision,  the  theoretical  result 

(42) 
obtained  by  Wadati  and  Toda  for  the  KdV  equation.      We  should  note  that 

the  above  mentioned  phase  shift  is  considered  in  (n,t)  space.   Since  the 

lattice  itself  is  contracted  toward  the  center,  the  phase  shift  is  less 

obvious  in  (x,t)  space.   The  phase  retardation  for  the  harmonic  lattice, 

(21) 
reported  by  JPW,  is  due  to  this  difference   '  ;  there  is  no  phase  shift  in 

(n,t)  space. 

The  phase  shift  observed  when  the  solitary  wave  hits  the  fixed  boundary 

(5.2.1)  can  be  explained  in  terms  of  the  two  wave  collision  process.   If  two 

solitary  waves  of  equal  amplitude  travel  in  the  opposite  direction  in  such 
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(a)      Head-on   collision. 
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(b)   Take-over  collision 


Fig.  5.13.   Collision  process  of  two  solitary  waves 
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a  way  that  a  collision  takes  place  exactly  at  the  location  of,  say,  m 
atom,  the  m   atom  never  moves.   Hence  the  fixed  boundary  condition  is 
strictly  satisfied  (Fig.  5.14a).   We  can  further  conclude  from  the  numeri- 
cal results  in  5.2.1  that  the  solitary  waves  lose  no  appreciable  energy  in 
a  collision  process  which  is  equivalent  to  the  fixed  boundary  condition. 
5.3.2   Energy  Loss 

We  can  extend  the  above  discussion  to  consider  a  different  initial 
phase  relation  of  two  waves.   Suppose  that  a  collision  takes  place  in  the 
middle  of  two  atoms.   Such  type  of  collision  is  equivalent  to  the  reflection 
of  the  solitary  wave  at  the  fixed  boundary,  except  that  the  interaction 
force  at  the  boundary  is  twice  as  strong  as  that  at  all  other  places  (Fig. 
5.14b).   Since  the  lattice  becomes  inhomogeneous ,  we  can  expect  that  the 
solitary  wave  may  be  deformed.   We  performed  a  numerical  experiment  on  the 
two  wave  collision  process  along  this  line,  and  discovered  that  some  of  the 
energy  is  actually  left  at  the  site  of  a  collision  (Fig.  5.15).   There  are 
50  atoms  in  the  lattice,  two  end  atoms  having  the  initial  velocities,  .5 
and  -.5.   A  collision  takes  place  between  n=25  and  n=26.   The  amount  of 
the  energy  taken  from  the  solitary  waves  is  very  small  (.37o  of  the  initial 
energy),  but  this  is  far  above  the  numerical  error.   This  numerical  experi- 
ment was  performed  twice  with  different  integration  steps  (i.e.  .02  and 
.01);  both  gave  the  same  results  within  the  numerical  errors.   Hence  the 
energy  loss  is  not  due  to  the  numerical  integration  scheme.   This  newly 
created  wave  consists  of  high  frequency  components  and  is  assumed  to  spread 
out  in  the  entire  lattice  at  t  =  °°.   Obviously  the  continuum  approximation 
cannot  treat  such  waves  properly. 

This  fact  may  seem  contradictory  to  the  theoretical  results  because  in 

(32) 
Toda ' s  lattice  the  only  effect  of  a  collision  is  a  phase  shift.       In 
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Fig.  5.14a.   Fixed  boundary 


Fig.  5.14b.   Another  type  of  collision  process 
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Fig.  5.15.   Energy  loss  in  a  collision  process 
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fact,  the  theory  states  that  two  solitary  waves  are  shown  to  exist  asymp- 
totically (t=+°°)  .   But  there  are  some  other  terms  in  the  expanded  form  of 
Toda's  solution,  which  vanish  asymptotically.   These  terms  are  neglected  in 
the  theoretical  argument.   Note  that  these  terms  arise  in  the  process  of 
approximating  the  solution,  not  the  equation.   In  our  numerical  experiments 
two  solitary  waves  are  created  at  both  ends  of  the  lattice,  between  which 
there  is  no  disturbance  initially.   Therefore  the  initial  conditions  are 
different  in  these  two  cases.   As  we  mentioned  before  the  wave  created  in 
a  collision  process  will  disappear  as  t  -*•  °°,  if  we  assume  the  lattice  of 
infinite  length;  asymptotically  there  are  only  two  solitary  waves. 

The  solitary  waves  in  all.  collision  experiments  propagated  after  the 
collision  without  any  decay.   No  decrease  of  the  propagation  velocity  due 
to  the  energy  loss  was  noticeable.   We  will  show  in  §5.6  that  the  energy 
loss  in  a  collision  is  more  obvious  for  the  hard  sphere  type  force  function 

In  the  thermal  situation  collisions  can  take  place  at  any  location  of 
the  lattice  with  equal  probability.   Also  the  energy  loss,  although  exists, 
is  negligible  in  the  small  sized  lattice.   Therefore  we  can  summarize  that 
the  only  important  effect  of  the  collision  in  our  model  lattice  (JPW  lat- 
tice) is  that  the  propagation  velocity  of  the  solitary  waves  is  effectively 
increased  due  to  the  phase  shift. 
5.3.3   Effect  of  the  Disturbances 

There  is  a  question  as  to  the  propagation  of  the  solitary  wave  when 
the  disturbances  (that  is,  the  incoherent  motion  of  the  lattice)  exist. 
Since  the  solitary  wave  is  in  the  sea  of  disturbances  under  such  circum- 
stances, it  is  very  difficult  to  isolate  the  energy  in  the  solitary  wave 
part  and  to  measure  its  change.  The  only  possible  way  is  to  observe  the 
propagation  velocity.   As  we  discussed  in  §5.1,  the  solitary  wave  is 
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uniquely  characterized  by  its  velocity.   Therefore  any  decrease  of  the 
velocity  means  that  the  energy  is  being  lost.   Fig.  5.16  shows  an  example. 

this  figure  was  obtained  by  JPW,  but  it  was  not  publicized  in  their 

*  (21) 

report  .       As  we  can  see  in  Fig.  5.16,  the  velocity  change  is  mainly  due 

to  collisions;  the  original  velocity  is  regained  after  each  collision. 
Hence  we  conclude  that  the  interaction  between  the  solitary  waves  and  the 
disturbances  is  negligible  for  the  lattices  of  order  100  atoms  or  so. 

This  does  not  mean,  however,  that  no  such  interaction  is  possible.   For 
example,  if  the  lattice  size  is  quite  large  the  solitary  wave  may  lose 
significant  amount  of  its  energy  due  to  the  interaction  with  the  distur- 
bances   as  well  as  due  to  many  collisions  with  other  solitary  waves.   As 
for  the  solitary  waves  of  very  small  amplitude,  we  do  not  have  any  way  to 
distinguish  them  from  the  disturbances. 

5.4   Time  Averaged  Kinetic  Energy 

We  define  this  quantity  as  follows: 

M 

K  (t)  -  7  >   y  2(knt)  (5.3) 

n       t  ^_ ,   n 

k=l 

where  M  is  the  number  of  time  steps  and  t  =  MAt.   This  quantity  gives  the 

local  temperature  when  the  lattice  interacts  with  the  heat  reservoirs  and 

M  ■*  ".   We  calculated  K  (t)  for  the  impulse  type  initial  condition.   We 

n 

considered  both  the  harmonic  and  the  anharmonic  lattices. 


The  author  is  indebted  to  E.  A.  Jackson  for  this  figure. 

**  r  1 

Saito  e_t_  al_.  reported  that  a  significant  decrease  of  the  velocity  was 
observed  for  the  solitary  wave  in  Toda ' s  lattice .  \^^> 
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5.4.1  Harmonic  Lattice  (Dispersive  Wave  Packet) 

As  we  discussed  in  Chapter  2  the  created  wave  packet  is  very  dis- 
persive.  That  is,  each  component  of  the  normal  modes  propagates  at  its 
own  group  velocity.   Consequently  K  (t)  is  a  decreasing  function  of  n  for 
a  fixed  time  t  (Fig.  5.17a). 

5.4.2  Anharmonic  Lattice  (Solitary  Wave  Plus  Tail) 

The  solitary  waye  gives  a  uniform  distribution  as  far  as  it  has 
propagated.   The  tail,  on  the  other  hand,  gives  a  decreasing  contribution 

to  K  (t)  due  to  its  dispersiveness .   The  resulting  distribution  is  shown 

n  ° 

in  Fig.  5.17b  (V.=.l)  and  in  Fig.  5.17c  (V.=.05).   The  latter  is  a  non- 
separating  wave  packet  case.   A  sharp  drop  of  K  (t)  in  Fig.  5.17b  shows  the 
location  of  the  solitary  wave. 

5.4.3  Temperature  Holes 

We  discovered  that  K  (t)  for  a  two  wave  collision  process  has  a  dip 

n  r  r 

around  the  site  of  the  collision  (Fig.  5.18a,b).   This  is  partly  due  to 
the  concentration  of  the  spring;  the  kinetic  energy  is  transfered  to  the 
potential  energy  when  two  waves  collide.   Obviously  this  effect  is  inde- 
pendent of  the  type  of  the  force  function,  and  in  fact  the  same  thing 
happens  in  the  harmonic  lattice  (Fig.  5.18c,d).   For  the  purpose  of  removing 
this  spurious  dip  from  K  (t)  we  experimentally  defined  the  temperature  in 
terms  of  the  total  energy.   Although  the  temperature  hole  disappears  in 
the  harmonic  lattice  with  the  new  definition  of  the  temperature,  it  still 
remains  in  the  anharmonic  lattice  (Fig.  5.18e,f).   Hence  we  must  conclude 
that  the  phase  shift  of  the  solitary  waves  is  also  responsible  for  the  tem- 
perature hole;  both  waves  are  effectively  accelerated  during  the  collision 
process,  contributing  less  to  the  local  temperature.   In  the  thermal  situa- 
tion the  temperature  holes  can  occur  at  any  location  of  the  lattice  with  an 
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equal  probability.   Therefore  it  will  not  leave  any  effect  on  the  tempera- 
ture distribution.   We  will  use  the  kinetic  temperature  for  the  actual  com- 
putation of  the  temperature  distribution. 

5.5  Piecewise  Linear  Force  Function 

As  we  discussed  in  Chapter  3  the  piecewise  linear  force  yields  the 

solitary  waves  in  the  continuum  approximation.   We  will  show  in  this  section 

that  such  solutions  also  exist  in  the  anharmonic  lattice  (in  the  discrete 
system).   We  choose  the  force  function  as: 

f(x)  =  x      x  >  -b 

(5.4) 
f(x)  =  2x+b    x  •  -b  . 

It  is  obvious  that  the  amplitude  must  exceed  the  threshold  value  b;  other- 
wise the  wave  motion  would  be  pure  harmonic. 

The  numerical  results  are  given  in  Table  5.2  for  various  b  and  the 
initial  velocities.   A  typical  profile  is  given  in  Fig.  5.19.   The  two 
wave  collision  process  was  also  examined  for  a  typical  value  of  V.  and  b. 
Two  solitary  waves  were  verified  to  exist  even  after  the  collision.   There 
was  no  phase  shift  nor  energy  left  at  the  site  of  the  collision. 

5.6  Hard  Sphere  Plus  Harmonic  Force  Function 

The  hard  sphere  plus  harmonic  force  function  first  introduced  by  North- 

.    (13) 
cote  and  Potts  may  be  considered  as  a  special  case  of  55.5.      This  model 

combines  the  harmonic  lattice  and  a  very  strong  nonlinearity ,  hence  it  is 

theoretically  very  convenient.   We  define  the  force  function  as: 

f  (x)  =  x       x  '  -b 

f(x)  «  -«      x    -b 
where  b  is  the  collision  distance  (Fig.  5.20).   Computationally  the  hard 
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Table  5.2 
Solitary  Waves  in  Piecewise  Linear  Force  Model 


b 

V. 

i 

Amp . 

V 

.01 

2 

.524 

1.35 

.05 

2 

.628 

1.33 

.1 

2 

.731 

1.29 

.2 

2 

.851 

1.25 

.4 

2 

.981 

1.16 

.6 

2 

1.025 

1.086 

.2 

1 

.491 

1.167 

b  :   break  point 

V.  :   initial  velocity 
l 

Amp.:   amplitude  of  the  solitary  waves 
V:   propagation  velocity 
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Fi  .  5.L9.   Solitary  wave  in  piecewise  linear  force  lattice 
f(x)  =  x       x  >  -  .2 
f(x)  =  2x+  .4   x  ■  -  .2 
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sphere  type  interaction  can  be  realized  by  interchanging  the  velocities  of 
two  atoms  when  a  collision  takes  place.   We  can  obtain,  on  one  hand,  the 
harmonic  lattice  by  increasing  b,  on  the  other  hand,  we  can  obtain  the  hard 
sphere  gas  model  by  letting  b  close  to  0.   In  our  numerical  experiments  we 
specifically  chose  b  =  .5.   Judging  from  the  results  in  the  previous  sec- 
tion we  can  expect  the  solitary  wave  solution  in  this  model  too.   The 
numerical  results  are  given  in  Table  5.3.   A  typical  profile  of  the  wave 
is  also  given  in  Fig.  5.21.   We  found  that  the  solitary  wave  is  quite  stable, 
provided  that  the  initial  impulse  is  sufficiently  large.   Note  that  the  wave 
motion  is  harmonic  if  the  amplitude  is  less  than  b.   For  some  critical 
values  of  amplitude  the  once  formed  solitary  waves  decay  as  they  travel 
along. 

We  also  examined  the  two  wave  collision  process  in  this  model.   Fig. 
5.22  shows  an  example.   In  this  example,  the  initial  velocities  are  .9  for 
both  ends.   The  following  facts  are  remarkable: 

(1)  the  energy  is  left  at  the  site  of  a  collision.   The  amount  of 
this  energy  is  considerably  large  (~4%) ; 

(2)  both  waves  are  advanced  after  the  collision; 

(3)  the  velocities  become  smaller  after  the  collision  due  to  the 
loss  of  the  energy. 

It  is  notable  that  the  energy  loss  in  the  collision  process  is  large 
enough  to  affect  the  propagation  velocity  in  this  model.   Presumably  the 
same  thing  could  happen  in  the  JPW  lattice;  the  nonlinear  effect  is  too 
small  to  be  noticeable.   It  is  quite  conceivable,  however,  that  the  soli- 
tary waves  encounter  many  collisions  in  the  thermal  situation,  if  the  lat- 
tice size  is  sufficiently  large.   Then  the  loss  of  the  energy  will  be 
significant . 


Table  5.3 
Solitary  Waves  in  Hard  Sphere  plus  Harmonic  Lattice 
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1.19 
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.5 
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1.0 
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.4043 
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.01502 

.02  5 
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.00248 
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Collision  of  the  two  solitary  waves  of  the  same  size. 
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Decay  at  n=38, 
b:   break  point 


V.:   initial  velocity  (initial  energy  =  .5-V.  ) 
i  ■  i 

V:   propagation  velocity 


E   ,  :   energy  in  the  solitary  wave  part 
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E   .  •   energy  in  the  tail 
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E   , :   energy  loss  at  a  collision 
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h:   integration  step 
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Fig.    5.20.      Hard    sphere   plus    harmonic    force, 
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Fig.  5.21.   Solitary  wave  in  hard  sphere  plus  harmonic  force  lattice, 
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Fig.  5.22.   Head  on  collision  for  hard  sphere  plus 
harmonic  force  model. 
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CHAPTER  6.   HEAT  CONDUCTION  PROBLEM  IN  ANHARMONIC  LATTICES 
AND  ITS  PHENOMENOLOGICAL  INTERPRETATION 


We  will  consider  the  heat  conduction  problem  in  this  chapter.   Our 
main  purpose  is  to  investigate  the  applicability  of  the  idea  of  the  soli- 
tary waves  to  this  problem.   As  we  discussed  in  the  previous  chapter  the 
JPW  lattice  and  the  hardsphere  plus  harmonic  lattice  (to  be  called  H+H 
lattice)  both  have  the  solitary  wave  solutions.   The  heat  conduction  problem 

has  been  numerically  considered  by  several  authors  and  the  temperature 

(21)  (44) 

gradient  has  been  observed  in  both  types  of  lattices.     ' 

We  can  obtain  a  better  insight  into  this  problem  by  locating  these  two 
types  of  lattices  between  two  extreme  cases;  the  harmonic  lattice,  the  JPW 
lattice,  the  H+H  lattice  and  the  hard  sphere  lattice.   They  are  put  in  the 
increasing  order  of  nonlinearity ,  in  a  sense.   The  last  model  will  need 
some  introduction.   In  this  model  each  atom  of  the  lattice  can  travel  freely 
in  the  interatomic  space  until  it  collides  with  its  neighbors.   The  colli- 
sion is  assumed  to  be  elastic  as  is  the  case  for  H+H  lattice.   Since  the 
velocities  of  two  atoms  are  exchanged  in  each  collision,  the  energy  can 
propagate  through  the  lattice  as  if  it  were  carried  by  one  atom  which  passes 
through  others  without  any  interaction.   Therefore,  if  this  lattice  inter- 
acts with  two  heat  reservoirs  of  different  temperatures,  the  internal 
temperature  distribution  will  be  uniform,  although  there  is  a  heat  flux  from 
the  high  temperature  side  to  the  low  temperature  side.   This  model  is 

equivalent  to  Knudsen  gas  model,  in  which  each  atom  literally  propagates 

(45) 
from  one  reservoir  to  another  without  collision.      The  harmonic  lattice, 

on  the  other  hand,  also  has  the  same  thermal  characteristics  as  the  hard 

sphere  lattice.   In  this  case  each  Fourier  component  (or  the  normal  mode) 

carries  the  energy  independently  of  other  modes  at  its  group  velocity,  hence 
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there  is  no  temperature  gradient.   We  also  derived  the  same  tendency  for 
the  semi-infinite  harmonic  lattice  (§2.6).   It  is  interesting  to  note  this 
conceptual  resemblance  in  these  two  models.   If  we  are  to  consider  the 
solitary  waves  as  the  media  which  carry  the  heat  in  the  JPW  lattice  and 
also  in  the  H+H  lattice,  we  encounter  a  serious  paradox:   the  solitary  waves 
are  conceptually  identical  to  the  individual  atoms  in  the  hard  sphere  lat- 
tice or  the  Fourier  components  in  the  harmonic  lattice.   Hence  there  would 
be  no  temperature  gradient.   This,  of  course,  is  a  contradiction  to  the 
observed  temperature  gradient  both  in  numerical  computation  and  in  the  real 
crystals!   We  will  propose  in  this  chapter  a  phenomenological  theory  to  solve 
this  apparent  paradox.   Our  discussion  will  be  mainly  qualitative  due  to  the 
following  drawback  in  the  present  stage: 

(1)  No  concrete  theory  has  been  established  for  the  solitary  waves 
in  the  JPW  lattice.  Therefore  we  must  use  the  empirical  rela- 
tions such  as  those  obtained  in  Chapter  5. 

(2)  The  stochastic  nature  of  the  present  problem  greatly  limits  the 
accuracy  of  the  solutions  due  to  the  fluctuation  therein.   A 
tremendous  amount  of  computation  time  is  required  to  obtain  the 
numerical  results  which  are  accurate  enough  for  the  quantitative 
discussions . 

In  the  first  part  of  this  chapter  we  will  analyze  the  numerical  results 
for  our  model  lattice  (N=50,  JPW  type  force  function).   In  §6.1  we  will  con- 
sider the  cooling  process.   The  heat  reservoirs  have  zero  temperature,  so 
that  they  work  as  heat  sinks.   We  will  compare  our  heat  reservoirs  and  the 
viscous  type  heat  sink  introduced  by  Nakazawa.      In  §6.2  we  will  adversely 
consider  the  heating  process  of  a  cold  lattice.   Only  the  anharmonic  lattice 
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will  be  treated  here  since  our  purpose  is  to  observe  the  formation  of  the 
temperature  gradient.   In  §6.3  we  will  consider  the  heat  conduction  process 
which  we  briefly  introduced  in  Chapter  1.   We  will  compare  different 
boundary  conditions.   In  §6.4  we  will  interpret  the  heat  conduction  phenom- 
ena in  terms  of  the  solitary  waves.   In  §6.5  we  will  discuss  H+H  model  and 
we  will  imply  the  importance  of  the  collision  process  with  regard  to  the 
formation  of  the  temperature  gradient. 

6.1  Cooling  Process  of  Anharmonic  Lattices 

We  computed  the  residual  energy  in  the  lattice  when  two  reservoirs 
have  zero  temperature  to  demonstrate  that  our  heat  reservoirs  actually  have 
viscous  effect.   We  also  examined  the  viscous  type  heat  reservoir  used  by 
Nakazawa  for  comparison  (§4.3).   In  the  latter  case  the  end  atoms  lose 
their  energy  no  matter  where  they  are  located  (§4.3).   The  results  are  in 
Fig.  6.1.   The  curve  I  corresponds  to  our  model  heat  reservoirs  and  the 
curve  II  corresponds  to  the  viscous  type.   The  initial  conditions  are.  the 
same  for  both  cases.   The  energy  loss  is  faster  in  the  latter  case  due  to 
the  continuous  interaction.   Our  model  may  be  described  as  having  the  infinite 
viscosity  inside  the  reservoirs  but  none  outside. 

6.2  Heating  Process  of  Anharmonic  Lattices 

The  purpose  of  this  numerical  experiment  is  to  observe  how  the  tempera- 
ture behaves  in  the  early  stage  of  the  heating.   The  only  difference  between 
this  experiment  and  those  in  the  next  section  is  the  initial  condition.   The 
temperature  distributions  are  given  in  Fig.  6.2.   The  reservoirs  have  the 
temperatures  .02  and  .01,  respectively.   At  t=0  two  solitary  waves  were 
produced.   Fig.  6.2a,  for  example,  is  the  temperature  distribution  at  t=50. 
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Fig.  6.1.   Cooling  process 

Curve  I  :   Impact  type  reservoirs 
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Fig.  6.2.   Temperature  distributions  in  warming  process 
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By  this  time  two  solitary  waves  have  interacted  with  the  reservoirs  and 
they  have  been  reflected.   Their  locations  are  at  n=8  and  n=35.   We  can 
clearly  see  the  individual  contributions  of  the  solitary  waves  and  the 
tails  to  the  temperature  distribution.   There  is  also  a  temperature  hole 
at  n=27,  indicating  a  collision  which  occurred  there.   As  time  goes  on, 
more  interaction  with  the  reservoirs  takes  place  and  the  temperature  rises. 
Sometimes  the  temperature  has  a  positive  gradient,  but  this  is  a  temporary 
effect  caused  by  the  large  amplitude  solitary  waves  which  happened  to  be 
near  the  low  temperature  end.   Such  an  irregularity  will  be  removed  with 
time.   The  effect  of  the  temperature  holes  also  disappears  with  time. 

The  numerical  experiment  in  this  section  is  preliminary  to  those  in 
the  next  section,  where  we  will  consider  the  long  time  behavior  of  the 
lattices . 

6.3   Heat  Flux  and  Temperature  Distribution  under  Different  Boundary 

Conditions 

In  this  section  we  consider  the  heat  conduction  problem.   We  performed 
numerical  experiments  with  different  boundary  conditions.   As  we  described 
in  Chapter  4,  the  location  of  the  end  atoms  were  checked  at  every  m  inte- 
gration steps,  where  the  integration  step  is  .05.   Specifically,  m=l,  5  and 
20.   Initially  the  atoms  were  given  random  velocities  sampled  from  the 
probability  density  (1.4)  with  T=.5(T  +T  ),  where  T   is  the  temperature  of 
the  left  reservoir  and  T~  is  the  temperature  of  the  right  reservoir. 
Specifically,  T  =  .02  and  T  =  .01.   We  ran  one  case  with  T  =  .012  and 
T„  =  .008.   The  latter  case  was  also  chosen  by  JPW.   In  some  cases  the 
atoms  also  had  the  initial  disturbances  in  the  displacement  too.   We  also 
ran  a  harmonic  case  to  single  out  the  qualitative  effects  of  nonlinearity . 
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Due  to  the  Limitation  of  the  available  computation  time,  the  maximum 
simulation  time  was  6000  dimensionless  seconds  and  in  most  cases  3000 
seconds.   Therefore  the  results  obtained  by  JPW  are  superior  to  ours.   The 
results  are  in  Table  6.1.   The  accumulated  flux  and  the  temperature  distri- 
butions are  also  given  in  Fig.  6.3  and  Fig.  6.4  for  two  typical  cases. 

We  will  summarize  the  important  facts  of  our  results  in  the  following; 

(1)  The  heat  flux  becomes  stationary  at  an  early  stage.   We  calculated  the 
heat  flux  from  the  data  on  the  accumulated  flux  by  the  least  square  fitting 
after  the  temperature  of  the  25   atom  becomes  steady.   In  the  long  run  the 
flux  thus  obtained  will  match  (accumulated  flux)/ (simulation  time).   The 
fact  that  the  flux  is  immune  to  the  temperature  fluctuation  (to  be  dis- 
cussed later)  seems  to  support  the  idea  of  the  solitary  waves  as  the  heat 
carriers,  since  the  solitary  waves  are  also  immune  to  the  background  dis- 
turbances.  We  can  consider  that  the  data  on  the  flux  is  quite  reliable. 

(2)  There  is  a  tendency  for  the  heat  flux  J  to  increase  with  m. 

(3)  The  temperature  gradient  is  a  very  fluctuating  quantity.   It  is  very 
difficult  to  obtain  reliable  data.   JPW  also  reported  this  difficulty. 
There  is  a  tendency,  however,  for  the  temperature  gradient  to  increase  with 
m. 

(4)  The  empirical  thermal  conductivity  k  =  J/(-7T)  is  dependent  on  the 
boundary  condition.   Since  the  fluctuation  of  the  temperature  gradient  is 
large,  we  cannot  make  quantitative  analysis  on  this. 

(5)  We  recorded  the  entire  history  of  the  interaction  between  the  lattice 
and  the  reservoirs;  that  is,  time  of  interaction,  incoming  velocities  and 
outgoing  velocities.   From  these  data  we  could  fairly  well  identify  the 
propagation  of  the  solitary  waves  of  large  amplitudes.   More  precisely,  for 
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Table  6.1 
Numerical  Experiments  on  Heat  Conduction  Problem 


Exp.  No. 

m 

time 

ci 

C 
r 

J(io"3) 

-7T(10"5) 

H 
T 

CJ 

T  /T 
1/2 

1 

1 

3000 

230 

267 

1.15 

1.56 

73.7 

1600 

.02/. 01 

2 

1 

6000 

221 

267 

1.06 

1.78 

59.6 

2000 

.02/. 01 

3 

5 

3000 

246 

325 

1.25 

4.4 

28.4 

1600 

.02/. 01 

4 

5 

3000 

245 

315 

1.23 

3.76 

32.7 

1600 

.02/. 01 

5 

5 

4000 

234 

294 

1.19 

3.94 

30.2 

2100 

.02/. 01 

6 

20 

2000 

262 

331 

2.00 

8.90 

22.5 

1000 

.02/. 01 

7 

1 

3000 

219 

269 

.434 

1.28 

34 

1600 

.012/. 008 

Note 


m:   boundary  crossing  of  the  end  atoms  is  checked  at  every  %   inte- 
gration steps. 

time:   maximum  simulation  time. 

C.  ,C  :   collision  frequencies  at  the  boundaries  (per  1000  sec). 

J:   averaged  heat  flux. 

VT:   temperature  gradient. 

h  •   empirical  thermal  conductivity  (-J/^T) . 

t  :   the  time  at  which  T   becomes  stationary. 

T  /T  :   the  temperature  of  the  reservoirs. 
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Fig.  6.3a.   Accumulated  flux  (Data  No.  2) 
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a  vel 


ocity  v1  given  to  the  atom  n=l  at  t=t   we  can  find  a  corresponding 


interaction  in  the  other  end  of  the  lattice  at  t=t  ,  where  t  -t   is  close 
to  the  time  it  takes  a  solitary  wave  to  cross  the  lattice  with  the  initial 
velocity  v„  •   We  used  the  initial  velocity  vs  propagation  velocity  rela- 
tion in  Fig.  5.6.   We  also  found  that  the  interaction  time  is  not  distributed 
uniformly  but  more  dense  when  a  solitary  waves  arrive  at  the  end  of  lattice. 
It  is  clear  that  the  solitary  waves  introduce  interaction  with  the  reservoirs. 

(6)  Histograms  of  the  velocities  of  the  end  atoms,  at  the  time  of  interac- 
tion with  the  reservoirs  were  obtained  from  the  above  data.   They  are  in 
Fig.  6.5.   The  distributions  after  the  interactions  are  merely  the  sampling 
from  the  probability  density  function  (1.4)  with  T=.02  for  n=l  and  T=.01  for 
n=50,  respectively.   We  could  not  determine,  however,  the  nature  of  the 
distributions  for  the  incoming  velocities  due  to  fluctuation.   This  diffi- 
culty was  also  reported  by  JPW . 

(7)  The  collision  frequency  is  higher  on  the  low  temperature  side.   This 
difference  is  necessary  to  retain  the  lattice  between  the  two  reservoirs. 

In  other  words  the  total  momentum  of  the  lattice  is  balanced  to  zero  in  the 

(21) 
time  average.   Consequently  the  pressure  is  the  same  on  both  sides. 

(8)  The  temperature  distribution  showed  a  periodic  structure.   This  phenom- 

(44) 
enon  was  also  observed  by  JPW  and  by  Helleman.       (Fig.  6.6.) 

(9)  The  heat  flux  J  is  considerably  lower  in  the  harmonic  lattice.   This 
result  was  also  reported  by  JPW  and  by  Helleman.   Since  the  collision  fre- 
quency is  also  low  in  the  harmonic  lattice  we  cannot  single  out  the  nonlinear 
effect  from  this  result.   We  can  explain  the  low  collision  frequency  in  the 
harmonic  lattice  as  follows;  firstly  the  anharmonic  lattice  is  stiffer  than 
the  harmonic  lattice,  therefore  the  "sloshing"  effect  of  the  lattice  is  less 
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Fig.  6.6.   Periodic  behavior  of  internal  temperature  distribution 
for  harmonic  lattice. 
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important  in  the  anharmonic  lattice.   Secondly  the  solitary  waves  propagate 
faster  than  the  sound  velocity,  resulting  in  more  collisions  per  unit  time. 
From  the  above  reasons  we  can  say  that  the  thermal  contact  is  better  in 
the  anharmonic  lattice.   In  the  present  scheme  of  the  interaction  between 
the  lattice  and  the  reservoirs  the  heat  flux  seems  to  be  limited  by  the 
boundary  effects  rather  than  by  the  properties  of  the  lattice  itself. 
Especially  the  "sloshing"  effect  of  the  lattice  is  not  likely  to  happen  in 
the  real  crystals. 

6.4   Phenomenological  Interpretation  of  Heat  Conduction  in  JPW  Lattice 

We  consider  the  heat  conduction  problem  in  the  JPW  lattice  in  terms  of 
the  solitary  waves  and  the  tails.   First  we  should  recall  that  the  inter- 
action between  the  lattice  and  the  reservoirs  produce  either  the  solitary 
wave  and  the  tail  (when  V.  is  sufficiently  large)  or  a  wave  packet  which 
does  not  yield  the  solitary  wave  and  the  tail  within  the  limit  of  the  lat- 
tice size  (when  V.  is  small)  (§5.1) .   The  numerical  results  in  the  previous 
section  as  well  as  in  Chapter  5  strongly  indicate  that  the  solitary  wave 
carries  the  major  portion  of  the  heat.  Because  of  the  high  concentration  of 
the  energy  and  the  high  propagation  velocity,  the  energy  is  transfered  very 
efficiently  (§5.1,  §5.2).   The  solitary  wave  is  not  absorbed  completely  in 
the  reservoirs  but  partially  reflected  (§5.2).   At  the  same  time,  new  soli- 
tary waves  are  created  due  to  the  interactions  induced  by  the  incoming 
solitary  wave.   Therefore  the  solitary  wave  also  works  as  the  "exciter". 
Because  of  the  high  propagation  velocity,  the  collision  frequency  between 
the  lattice  and  the  reservoirs  is  considerably  increased  ((9),  56.3).   An 
important  fact  is  that  the  momentum  balance  of  the  entire  lattice  ((7), 
§6.3)  is  maintained  chiefly  by  the  "catch  ball"  of  the  solitary  waves  between 
the  two  reservoirs. 
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The  tail  (including  the  wave  packets  for  the  low  initial  velocities) 
is  characterized  by  dispersiveness ,  low  energy  and  low  propagation  velocity. 
Its  contribution  to  the  heat  flux  is  not  significant,  but  the  individual 
tail  gives  a  decreasing  contribution  to  the  temperature  in  the  direction  in 
which  it  propagates  (§5.4).   We  should  note  that  there  is  a  big  difference 
between  the  waves  in  the  harmonic  lattice  and  the  tails  in  the  anharmonic 
lattice,  although  both  waves  are  harmonic.   In  the  former  case,  the  harmonic 
waves  are  responsible  for  the  momentum  balance  of  the  lattice  as  well  as 
for  the  heat  flux,  whereas  the  tails  are  responsible  for  none  of  them  in 
the  latter  case.   In  other  words,  there  is  a  constant  (small)  flow  of  the 
tails  from  each  reservoir  in  the  anharmonic  lattice,  the  flow  rate  being 
determined  by  the  traffic  of  the  solitary  waves.   We  can  expect  that  the 
unbalanced  flows  in  two  directions  due  to  the  temperature  difference  will 
result  in  some  gradient  from  the  high  temperature  to  the  low  temperature 
side .   This  is  what  we  propose  to  solve  the  paradox  which  we  stated  in  the 
beginning  of  this  chapter,  since  both  the  collision  process  and  the  inter- 
action between  the  solitary  waves  and  the  disturbances  do  not  seem  to  be 
significant  in  our  model  lattice.   Unfortunately  we  do  not  have  enough 
numerical  evidence  as  yet  to  support  this  proposition  quantitatively. 

The  disturbances  in  the  lattice  are  the  collection  of  dispersed  tails. 
Since  the  tail  loses  its  coherency  as  it  propagates,  the  disturbances  are 
quite  incoherent  and  do  not  contribute  to  the  heat  flux.   The  energy  of  the 
disturbances  does  not  increase  indefinitely  but  is  maintained  at  a  certain 
level  on  the  average  due  to  the  interaction  with  the  reservoirs. 

The  heat  flux  and  the  temperature  gradient  are  related  to  each  other 
through  the  nonlinear  coefficients;  the  propagation  velocity,  fraction  of 
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the  energy  in  the  solitary  waves  and  the  stiffness  of  the  lattice  are  all 
dependent  on  the  nonlinear  coefficients. 

The  effects  of  the  boundary  conditions  on  such  quantities  as  the  heat 
flux,  the  temperature  gradient  and  the  thermal  conductivity  are  not  clear. 
Presumably  the  energy  is  more  efficiently  exchanged  at  the  reservoirs  if 
the  end  atoms  have  more  freedom  (i.e.  m  is  large).   Therefore  we  can  expect 
higher  flux  in  such  cases. 

At  any  rate  the  dependency  of  the  thermal  conductivity  on  the  boundary 
conditions  implies  that  this  quantity  is  not  intrinsic  to  the  lattice. 

Finally,  we  will  summarize  the  results  of  this  section. 

(1)  The  solitary  waves  carry  the  major  portion  of  the  heat. 

(2)  The  solitary  waves  considerably  increase  the  collision  frequency  with 
the  reservoirs,  due  to  high  propagation  velocity. 

(3)  The  solitary  waves  also  maintain  the  momentum  balance  of  the  lattice. 

(4)  The  tails  are  one  possible  source  of  the  temperature  gradient  in  the 
JPW  lattice. 

For  other  possibilities  of  the  source  of  the  temperature  gradient, 
refer  to  the  next  section. 

6.5  Hard  Sphere  Plus  Harmonic  Type  Model  and  Other  Possibilities  of  the 

Temperature  Gradient 

As  we  discussed  in  the  previous  section,  the  JPW  lattice  has  the 
characteristics  of  Knudsen  type  model.   We  next  consider  the  H+H  type  model, 
Helleman  investigated  the  heat  conduction  problem  in  this  model  and  dis- 
covered that  the  temperature  gradient  becomes  more  obvious  than  in  JPW 
lattice.      We  showed  in  §5.6  that  H+H  type  lattice  also  has  the  solitary 
wave  solution  and  we  discovered  that  a  non-negligible  amount  of  the  energy 
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was  lost  from  the  solitary  waves  in  the  two  wave  collision.   From  these 
results  we  can  consider  the  collision  process  as  very  important  in  the 
strongly  anharmonic  lattices.   Namely,  the  solitary  waves  lose  significant 
amount  of  energy  as  they  collide  with  other  solitary  waves.   In  this 
picture  the  relation  between  the  heat  flux  and  the  temperature  gradient 
becomes  inherent  to  the  lattice.   Another  possibility  is  the  interaction 
between  the  solitary  waves  and  the  disturbances.   In  the  JPW  lattice  the 
effect  of  the  disturbances  is  not  important.   But  if  we  increase  the 
number  of  the  atoms  enormously,  such  effect  may  not  be  negligible  any  more. 
In  this  case,  too,  the  solitary  waves  will  decay  as  they  propagate  through 
the  lattice.   No  research  has  been  done  in  this  area,  however.   We  will 
need  a  very  high  speed  computer  of  the  Illiac  IV  class  to  perform  numerical 
experiments  with  super  large-sized  lattice. 
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CHAPTER  7.   CONCLUSION 

In  this  thesis  we  considered  the  properties  of  the  solitary  waves  both 
analytically  and  computationally.   In  the  first  approach  we  introduced  the 
piecewise  linear  function  model,  by  which  most  of  the  special  characteris- 
tics of  the  solitary  wave  solution  and  also  the  periodic  solutions  can  be 
explained  with  the  elementary  functions.   As  for  the  second  approach,  we 
performed  numerical  experiments  on  the  anharmonic  lattice  introduced  by 
Jackson,  Pasta  and  Waters.   We  discovered  that  the  impulse  applied  to  the 
end  atom  of  the  lattice  produces  the  solitary  wave  very  efficiently.   By 
this  method  we  performed  various  kinds  of  experiments  as  to  the  properties 
of  the  solitary  waves.   One  remarkable  thing  is  that  the  solitary  waves 
thus  created  lose  their  energy  in  the  collision  process.   This  has  never 
been  reported.   Although  the  amount  of  loss  at  each  collision  is  very  small, 
this  fact  indicates  the  picture  of  the  decaying  solitary  waves  in  the  an- 
harmonic crystals.   We  could  also  show  that  the  piecewise  linear  function 
model  and  the  hard  sphere  plus  harmonic  type  model  also  have  the  solitary 
wave  solutions. 

Our  numerical  results  were  applied  to  interpret  the  heat  conduction 
phenomena.   Numerically  we  could  show  that  the  major  part  of  the  heat 
flux  is  carried  by  the  solitary  waves.   The  origin  of  the  temperature 
gradient  is  still  unclear,  but  the  dispersive  wave  packets  which  are  created 
with  the  solitary  waves  while  the  lattice  interacts  with  the  heat  reservoirs 
might  be  a  possible  source.   The  collision  process  and  the  interaction  with 
the  disturbances  seem  to  be  negligible  in  our  model  lattice.   Our  numerical 
results  indicate  that  the  heat  flow  in  our  model  lattice  is  near  Knudsen 
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type.   It  is  conjectured  that  the  thermal  conductivity  will  become  the 
intrinsic  quantity  by  increasing  the  lattice  size  considerably. 

Finally,  we  mention  that  the  numerical  integration  method  due  to 
Nystrom  was  found  to  be  superior  to  the  classical  Runge-Kutta  type  inte- 
gration method. 
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APPENDIX  A 
SOME  PROPERTIES  OF  TCHEHYCHEFF' S  ORTHOGONAL  FUNCTIONS 


Tchehycheff  s  equation 


,1  2^    d^      dy     2 

(1-w  )2"Wdt+ny  =  0  (A<1) 

dw 


.   .    (25),  (26) 
has  two  solutions 


and 


y   =  T  (w)  =  cos(narccos  w)      first  kind  (A. 2) 


y9  =  U  (w)  =  sin(narccos  w)      second  kind  .         (A. 3) 


Note  that  the  above  solutions  differ  from  the  usual  solutions  by  the 

factor  2     (according  to  van  der  Pol     ).   T  (w)  is  a  polynomial  of 

n   order,  whereas  U  (w)  is  the  product  of  Vl-w   and  a  polynomial  of  n-1 


order . 

Explicit  solutions  are  given  as  follows  for  n  =  0  ~  4 : 


TQ(w)  =  1  UQ(w)  =  0 


T  (w)  =  w  U  (w)  =  Vl-w 

T  (w)  =  2w2-l  U2(w)  =  2w7l-w' 


3  2    7   2 

T  (w)  =  4w  -3w  U3(w)  -  (4w  -l)v/l-w 

L.         9  3      I       2 

T  (w)  =  8w  -8w  +1  U4(w)  =  (8w  -4w)Vl-w   . 

Note  that  T  (w)  is  an  even  function  of  w  when  n  is  even,  whereas  U  (w) 
n  n 

is  an  odd  function  oi  w  when  n  is  even. 

TT 

Since   arcsin  w   =  —  -   arccos   w, 

cos[  (2n+l)arcsin   w]   =   cosf.  (n+2)n-  (2n+l  )arccos   w] 

=    (-1)    sin[ (2n+l)arccos   w] 

"    ('1)nU2n+l(w)    '  (A'4) 
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The  relation  between  Bessel  functions  and  Tchehychef f ' s  polynomials 

is  best  explained  by  using  Hansen's  integral  representation  for  the  Bessel 

(25) 
function  of  integer  order     : 

•    "n  n  -4-  ft 

,        X  1  f-ltCOSO  Ann  •  .  r-  N 

J  (t)  =  ——  J   e       cosn9  d9  .  (A. 5) 

o 

By  introducing  a  new  variable  w  =  cos9, 

.  -n   -1  .. 

/x    i     f     itw    r      -1  -i ,       .  Ax-1  , 

J  (t)  =  e    cosLncos   wj(-sino)   dw 

n       TT  « 

.-n   1   T  (w)   . 

1    ?  n     iwt,  ..  ,. 
I     ■     e    dw                        (A.  6) 

!    -1/7-2 

1-w 

The  above  relation  shows  that  the  Bessel  function  has  no  higher  frequency 

2 
components  than  w   >  1.   This  corresponds  to  the  fact  that  we  are  dealing 

with  the  low  pass  filter  (§2.1). 

Furthermore , 

-JV~  =  ^(J  ,l(t)+J   l(t)) 

t      zn   n+1      n-i 

. -n+1   1  iwt 

=  —: f    (-T   ,(w)+T   .  ) )   6     dw 

2nn  J    v  n+lv  '  n-ly/  / ~ 

•J  1-v 


-w 


. -n+1   1 

1      0  iwt 

I    U  (w)  e    dw  , 

nn  *J     nv 


since 


T   . (w)  -  T    (w)  =  cos[ (n-l)arccos6]-cos[ (n+l)arccos9] 
n-1       n+1 


=  2sin (narccos9)sin(arccos9) 


=  2U  (w)-Un  (w) 
n      L 


-ifi    ' 


w    U  (w)  . 
n 
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Hence, 


J    (t)  2.l-n      1 

— =  z~  U    (w)cos(wt)dw 

t  nn     <J         n 


(A. 7) 


when  n  is  an  odd  integer, 
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APPENDIX  C 

d   2 
ON  THE  TWO  SOLUTIONS  OF  (—)   =  P(w) 

The  differential  equations  which  we  are  interested  in  is  of  the 

following  form: 

,   2 
O   -  P(w)  (C.l) 

where  P(w)  is  a  cubic  or  quartic  polynomial. 

We  showed  in  Appendix  B  that  the  solutions  of  the  above  equation  are 

expressed  in  terms  of  the  Jacobian  elliptic  functions,  sn,  en  and  dn. 

Corresponding  to  the  fact  that  the  circular  functions  sine  and  cosine  are 

both  the  solutions  of 

,dw        2 
(dz}   =  1_W 

with  different  initial  conditions  at  z=0,  there  are  also  two  solutions  to 
(C.l)  with  different  initial  conditions  at  z  =  0.   These  two  solutions  can  be 
matched  by  shifting  one  of  them  by  a  quarter  period,  as  is  the  case  for 
the  circular  functions: 

TT 

cos  z  =  sin(z+  — )  . 

For  example, 

,   2 

(j1)   =  (a-w)(w-b)(w-c) 
dz 

has  two  solutions, 


/   u  \   2  ,va-c    ,  .        .2   a-b 

w  =  a-(a-b)sn  ( — —  z  k)   ,      k  =  

f  2  a-c 


and 


,   a-b   2  ,\a-e     . 

b  -  c s  n  (-* — r—  z  ,  k ) 

a-c       L 
w2  =  — 

l-(— )sn  (*-y  z,k) 
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By  using  the  following  formula 


(47) 


sn2  (x+K)  =  ^Sl   =  lZ§fJf 

dn  (x)    1-k  sn  (x) 


we  can  show  that 


WL  (z+  ~ <k)  =  a-(a-b) 

Va-c 


1-sn 


,   a-b   2 

1-  sn 

a-c 


a  (a-b)   2  2 

a-  — * sn   -  a+b+(a-b)sn 

a-c 

a-b   2 

1-  sn 

a-c 


,   (a-b)(a-a+c)   2 

b-  J " Lsn 

a-c 


,   a-b   2 

1-  sn 

a-c 


a-b   2 

b-c-a^c~Sn 

,   a-b   2 

1-  sn 

a-c 


w2(z,k)  . 


Hence  we  can  choose  one  of  the  solutions  to  (C.l)  to  discuss  the  nature 
of  the  periodic  solutions  without  losing  generality. 
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APPENDIX  D 
COMPUTER   PROGRAM 


r i ^e i   mil/xmv 
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tjF.  CjIN 
t%    THIS    PROGRAM    COMPUTES     THE.    HEAT    FLU*    ANT     TEMPERATURF     GRADIENT 
*%    OF    PME    UlMENSIUNAL    AmhARMIJNIC    LATTICE    wTTH    FORCE    FUNCTION 
t%  Fs**ALPHAX**2    ♦    RI*X**3 

X(AND    IMPACT     TYPE    HEAT    RFSERvOIRS    AT    bOTH    ENDS.* 
%%    wrITTEn    2/24/73    Ano    rEvIsl)    3/25/73         k«mI|JrA 

FILE     STATK)N(KlNF)st3#MA<RECSlZE«20,lNTMnDF»4.MYUSF»3)  J 

file  line(Kind»7»maxrecsiZe*22)j 
fURmat  nn T 

FLXXlO»F*.2»X3'ElO«3»X3»E10»3)» 
FHF(X6WV.2»XliElOt3tX3»E10.3)f 

FlKMt  13.  10(X?,R9.2))  i 
LABEL  L4.L5j 

REAL    ALPHU'ALl3»RI»Rl4»R2,Ht>'HZSCL'H2,H4»HIVHl5»ZMXJ 
REAL    M»*5. Kft, K7. K«, K9.H.  I3»  16.  19.  I192.HI19?»hM23»     H«o* 

HEAL    AMP,wnTH,SUM,TMP,RFLUX.LFLllX,TMH.YiPl,Yli,Zjl,H?i3» 

REAL  ARHA*    Y,Z[0»  5  1]  .  TP  .  T<5  [  1  «  50  .  0  «  5l]» 

REAL  ARHAY    K  1  , «2 , K 3[  1 : 5 1 ] . K NTp[ 1  I  50) , LF. RFT * » ^00 1 » 

INTEGER  AKRAY  XC  t) »  YCflt  1  «  1  00  JJ 

INTEGER    LCNTR.RCNTR.TW035' 

INTEGER     I.J,K,L.H,JJ.KK,LL»L1.FLG,FLG?.FLG3, 

INTEGER     X;      XFOR  RANDOM  NUMRERS 

INTEGER    M.TIME.STP.ISTAHT»ILAST»HI> 


HEAL  P 

H 

Xt=SECO 

(271 


REAL  P 

Reg 


ENn 

PRUCED 
VAL 

INT 
REA 
REA 

%%%%t%%%%ii^%%% 
RE 


RJCE 

EGIN 
NUrtO 
U26l 

Ni)» 
*i  ■** 

ROCE 

IN 

*UNE 

* 

REAL 

TMP« 

MaXW 

;  % 

Hit, 

URE 

UE 

EGER 
L  aR 
L 
%%%* 

GIn 
REAL 
INTE 
REAL 


i%t%%t%%%%%%ttt%ti*%l%%%%ii*l%%%%%t %%%%%%%%%%% 

UiIRE  RNDXj 
l%ti%%i.lZ%%%tt%l*%i%%%%%%%%l%%%%t%%%t%%%%1%%Xt 

<*  OUE  TU  D.KNUTH 
R0(0lNTEGER(3iai592653  MUX  X)+OTNTFGfR 
«29))  Mni)  Tw035» 

RNUX:»X/T^035; 
*RNnX  X* 
*%%*%% tt%%l%%it%t%*%%%%%%t%t%l%l%%%%X%t%%%%%ti 

DURE    MaXWelL(T)'   REAL   T» 

l%%t%%t% *%*%%%*%%%*% 1*%%%%%%%%%%%%%%%%%%%%%%%% 

OIMFNSIONAL  MAXWELL-HOLTZMANN  OlSTRlROTlON 

TMp: 
=-T*LN( 1-RNDX) » 

ELL,=S«RT(TMP+TMP);  %%%j     IS    ALWAYS    ASSUMfD    T 

*    mAxwEll"ROlTZmAnn    i% 
i*,t>t%i%%%*t%%%*%%%*%%%% 
LSS^fA.IS.lL»INTL»lCPT»SLHP); 

IS.IL.INTL; 

IS,  IL, INTLj 
RAV  A[* J  ; 

I C  P T  » S LO  P  J 
i%)k*%%%%%%l%l%%%lit%%%l 

Tl  J 
l»  E  H  I  .  J  I 

SUMY, SUMZ. MAX. TEMP; 
SUMYt*SUMZ:»MAX:sTEMP|sOJ 


OOOO0100 
00000200 
00000300 

oooooaoo 

00000500 
00000600 
00000700 

noooottoo 

00000900 
OOOO1O00 

OOOO1100 
00001200 

00001300 
OOOOHOO 

OOOOl^OO 
00001600 
00001700 
00001800 
00001*00 
00002000 
00002100 

OOOO2200 
00002300 
00002400 
00002500 
00002600 
00002700 
OOOO2800 
S***X**OOOO2900 
00003000 

***s***oooo3ioo 

00003200 
00003300 

oooo3aoo 

OOOO3500 

00003600 
*«**»**noOO3700 

OOOO3800 
*XXXX**OOOO3V00 

OOOOUOOO 

00004100 

00004200 
OOOO4300 
00004400 
0  RE  POOOOO4500 
00004600 

OOOO4700 
00004800 
00004900 
00005000 
00005100 
00005200 
00005300 
00005«00 
OOOO5500 
00005600 
OOOO5700 
00005800 
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fur  i»  =  is  step  intl  until  il  do 
begin 

J'=I/INTL' 

SUMYl.SUMy+I; 
TEMPJ=TEMP*1*A[ j  j i 
SUMZ««SUHZ*Al Jl I 
MAXt=MAX*I*I» 
ENDj 
j:=(  I|  -  I S)/ I NTL+ 1  I 
II        l=J*MAX"SUMY**2j  , 

Kpt:s(sumZ*max-sumy*temp)/ti* 


S*«lNlTERCFPT 


pRncLnuw 

VALUE 

JNTE^FR 

INTEGER 

^eg  In 

FORMAT  I) 

F2(X«),.f- 

F3<x5.«0 
ri<  (  X  3  »  1  3 

fl  W  N  I  N  T  E 
pnlNTFP 
VALUE  AH 
(4*04010 
'4 "  2  6  1  8  1  V 

»l   |<3>?(# 

n  "  a  » f  9  8  t 

4"A6A7AH 

K 

« 

PI  «PlJ  I  NT 

Ql=pUINT 

REPLACE 

WRITL(ST 
IF  SCAl  E 
RE'ilN 

SCALE 

WRITE(ST 

PflR    I  I «0 

IF    I 


EM„ 

F     HAZE 
N.MS; 
fJ.MS; 
AHKAY 


SLnP,*(J*TE«P-SUMY*i)UMZ)/Tl|        *»SL()PE 
Up     L^AST     SijyARE     pITTINt,; 

l.PLUTCXCIl,YCC3»N#MS)» 


XCi),  YLf)[  *  i  t 


UT     F0U"A1  1UA1  1  |)Al  1FA1  1  1040200")     . 
»    >    »U'    ••»•*;•     •  •  #  t  <    » 4    0  0    >  • 

fi ,  x  6  4  ,  n  -  n  ,aMooM). 
ii,lJ(i«i'f«j"),a».()0,»). 

(  j3»X2)'a"Ali200,'>» 

GtH   SCALE*- Lf,J  INTEGER  MARRAY   X  r  o  «  0  ]  » 

p.  q; 

RaY   ASCII 
2u3J72n2E2F16US?SOROCOf)OEOFlOlll?133C3o32 

3F2»'lClL)lElf"**  I  """*J»8'  (  )*♦•"./ 0  12345678 

AriCUEKr,HljKLMN(lpQRSTlJVwXYZ[,,a,'EO"«)-,'' 

u  .a   i<4..ai.M<.H/UMMoa.o..u->u..  u.  OiUrouou.. . 


Hi.dj(l(48b06rt<rM«tt99i9<!93V4Vt,96yr'9ttPV 
A9CUrSAD0AlU("'). 


9" 

A2A3A4A5" 


LLS 
WRITL(S 
WRITE(  S 
WRITL(S 
FNU, 

writer 

*RITE(S 

» 

FOR  II" 

REGIN 


E*(*)*?f 
ER(AScU)J 

PUINTF.Rl  Xl  U  J  )  HY  a"  A  1  1  1  ••  FUR  2  »  "  "  FUR  2» 

4 "At  no"  FUH  2' 
A  T  I  UN.FO  )  I 
F  L  lf  »  0  THEN 

FLfil«l  I 

A T I UN »  <« " A U C A 11 C  A 1 1 9  A 1 1 9 00">  i ' 

STEP     1     UNTIL    2A     OU    HEtilN 
m  U  L)    b  =  U     T  H  L  -j 

tfRlTE<STATIl.N*Fl»<15aI>/H£SCL> 

E        WHI TE(ST AT  IUN.F2)  i     END! 
TATIUN, f3  )  > 

tat  I  un  •  r  a  »FUN    l»«o   step   5   until   60   0 11    in 

TATIlJN.<4"AilFAllFU0">>l 

TATIUN.  <4"Al  1 200n>  )  I 

TATIUN.<nTlML«".Te>.X3."ENRc.Y*H.R11.5.X3»,,ALPHA«,'.En.4> 

MS»SuM. ALPHD ) J 

1  STEP  1  UN  UL  n  no 

YCUl  I  )  t*YCUl  1  )  +  \S  I 


REPLACE  P  HY  u*XCurIJ  EUR  l.u+YCOII]  EUR  1J 


00005900 
00006000 
00006100 
00006200 
00006300 
0000*400 
00006500 
00006600 
00006700 
0000*tt00 
00006900 
00007000 
00007100 

S**X*X**00007200 
00007300 
00007400 
00007500 
00007600 

*  %  *  %  %  I  %  t  0  0  0  0  7  7  0  0 
00007800 

0000/900 
00008000 
00008100 
00008200 

0000*300 
00008400 
00008500 
00008600 
00008700 
00008HOO 
O0O0A900 

00009000 
00009100 
00009200 
00009300 
00009400 
00009500 

0()009ftOO 
00009700 
00009800 
00009900 
0  0  010  0  00 
00010100 
00010200 
0  0  010  3  0  0 
0  0  010  4  00 

0  0  010  50  0 

0  001 06  00 
00010700 
00010800 
00010900 
0001 1000 

0001 1 100 

0001 1200 
0001 1 300 

0001  lion 

00011500 

0001  uoo 

0001 1 700 
0001 1800 
0001 1900 
00012000 
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WRITHSTATIUN*  liX[*l)l  0001210C 

end»  00012200 

write* stat  ion. <4"a112a11900">>»  000 12  300 

END    i)F    Ha-?ELTINE;  00012^00 

PROCEDURE     KINETICTEMP(N»X,TIME.N1  >'  00012600 

value  NtTiME.Ni>  00012700 

i  NTEtiLK       N.TIME.N1I  0001280C 

HEAL    AHRAY       *[*]>  0001290C 

*I«I*I«J»X*X%*«*******»«X**»**«»*4«»*4»*i»I«»X^8*I«Hil-«I«I«««I«««*«lO00  13000 

fit&TN 

FORMAT  out 
F51       1X5, w  KINETIC  TEHPERATORE  AT". 
"  TI*E  «".14)» 
E^2(  x20»"  I  «.  5(9(".")»"  |  ").  )  » 
E56(xl'El3.6»X2M2»X2'"0rt')» 
F54(Xl.El3.6.X2.l2.X2,"o".*<''X''),), 
FS7(X20»"0"»Xtt»".<!M»x8.,,.4,,»x8.,,.6,,.)<8,",8"»)<7. 

"1 ,0") | 
INTEGER    Itj.KI 
WEAL    T1'T2»  TAX'TEMPI 

REAL   ARRAY   FOUR  I  EH [ 1  ' N ] » 

W«ITE(LINE[SkIP  1])» 
TMH!=TIME*HH 
T  A  X  ;  =  A  M  P  ♦  A  M  P  > 
FUR  I J31  STEP  1  UNTIL  N  00 

TAXiBrtAxCCFUURlEKtl] l"X(I ]/TMH)*T*X)) 
wRITE(LINE    »f  51  . TIME  )  * 
WRITE(LINE    .F57); 
*RITE(LINE    »E52); 
Tlj=100/TAx; 
FUR  !»  =  *  STEP  *  UNTIL  N*       nO 
BEGIN 

K:ST1*(TEMP:3F0URIER[I1 )l 
If     a    NEO  0  THEN 

WRITE(LINE    »F5a»TEMP»I.K5  El SE 
WRITF(LINt    .F56.TLMP. I )> 
ENO, 
*<RITE(LINE    .F52); 
"R1TL(LINE    'fK>7); 
rtRITE(LINE.</«ACCUM.LFLUX  =  ".Rl?.5,xlJ,"ArcUM.  RfLUXs" 
»Kl2.t)/>'LFLUX»KFLUX)  J 
l.SSQ( FOURIER* 1  .  N  .  1  »  T  1  .  T  2  )  > 

*RlTE(LlNE,<"SLUPL=".E12.5.X3."LEFT  LSO  TMP»". 
E12.5,X3,nRlbHT  LSQ  TmPc".E1?,5>.t2,T1+T?, 

Tl*N*T2)l 

WRITE(LINE»</X10»"LEFT  COL="» I 6. X35. "R I GhT  cnL="»I6>. 
LCnTR.RcNTR); 
wrITE(  STATION. Fi»l»  TIME); 
wRITE( STATION, <"ACCUM.LFLUX=M.R 12. 5. x5."ArCUM.  RFLUXa" 
,R12.5/.MLEFT  COLs".  I6.X5."Rlr,MT  CHL  =  "  •  I  «  •  //  » 
"SLOPE    n.E12.5./»"T25    « .  E  1  2  .  5>  ,  LFL'Jx*  .  5  • 


00013100 

00013200 
00013300 
OOO13A0C 
OOO1350C 
OOO1360C 
OOO137O0 
00013800 
0001390C 
00014000 

oooiaioo 

OOOH200 
OOO143O0 

OOOlA^OO 
OOOH500 
0001460C 
OOO1470G 
0001480G 
00014900 
00015000 
00015100 
00015200 

00015300 
0001540C 
0001550C 

0001560t 
O001570C 
O00l580fc 
0001S90T 
000 16000 
00016100 

00016200 
00016300 
00016400 
00016500 

00016600 
00016700 
OO016800 
00016900 

00017000 
00017100 
OOO17200 
0001^210 
OOO17300 

oooi7aOO 


RFLUX*.5,LCNTH.RCNTR.T2, FOURIER [25])  I 

END;   «*XKINETIC  TEMPE-RATUKL    %%%%% 

%%%t%%%%ti,%%%%ti1,ii,%i%%%t%tii,%tn%%t,i,%i,%%%%%%i%%%t%%%%%%l%%%%t>%%%%%%%%%noOl750Q 

PROCEDURE   PnLYGRAPH(P. N.w. SCALE) >  00017600 

VALUE    P.N, SCALE;  OOO17700 

INTEGER    P,N. SCALE;  OOO178O0 

REAL  ARRAY  W[*.*l;  00017900 

%%%*%%%% t *%%%%%% it  tt%t.%%t%%%%%%%ttt%i%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%t  000  I  SOOO 
CqMMEnT      THIS  PRqCEDURE  OIsPLAYs  10  GRaPHsUsInG  A  LINE    00018100 
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begin 

FflHM, 


printed. they  shuuld  not  overlap. w  contains  the 

OATA.N  IS  THE  NUMtiEH  OF  GRAPHS.    P  tS  THE 
STARTING  POINT  UF  ARSCISSA.j 


F«( 


In 

In 

wKITF.d 


AT  OUT  FO(  x4.  12(  13. x7)  )  . 

»«  F2(  Xl  »  l4»Xl  .  W(  x*  ."*"))  i 

FJ(X6»il(,,X,'»*("."))»,'X,,)» 
FFF  (  I3.10CF12.8). I«)» 
[ii lU(Fli.lO)* I4)« 

F2(Xl.I4,Xl»X*»".".2(X*.n*"»X*."«i".X*. *•♦''. X*» 

"()"»X*'"»"))J 
TEtiER   t»J#K.TMP»L' 
rt*»t"R  ARRAY    VdUOJi 
lNt,<«vl».,xl2»"v2  ".xll."OUAuM»Xin,«CUnE»,X'n»  "SCALE 

xio,"stfpsUe">)  ; 

"HlTE(LlNL.<6(Rll.a#X3)/>,AMP.rt0TH.ALPH0.HI,SCALF»H)! 
FUR  l«sO  STEP   1  UNTIL   N  0 1  V  10  -1  DD 
BEGIN 

«RITE(LINE»FU'Ff)R  JJ=0  STFP  1  IJNTTL  10  DO 

►*♦  I  *  1  0*  J  J  » 

*RITE(LINE,F1 ); 
F  OR    J:sl     STEp     1     ijNTI  L    N  D" 

BEGIN 

FOR     K  l  « 1     STEP     1     UNTIL     10    m 

IF     TMPlsIMO  +  K    LE«    N    TMFN 
Yt K  ]  :  =    w[ TMp.  J)*SCALF» 
VRITE(LINL'F?»J  .U'^Ytll-'  +  Yt^TYn]' 

9*Y13 J-Y12 J»«+Yt4]-Yl1]»9+Y[5]-Yt4l« 

Vvyt  9]-Yl8J  .  V*y[  I  (J  ]  "Y  T  V  1  )  S 

END, 

^RITE(LINE.FJ) J 

«RITE(LINE  tSMP     ll.FO'FlIR    JtxO    STFP 
UII     (P*        J  *10*.l)  )  J 
FUB      u»0    STEP    50    UNTIL    N    0  I  V    50-1     On. 
MEr,  IN 

wRITL(LINE»<X3»10(I'3.X7)/>»FnR  L  I  =  1 

UNTIL  10  D0(P*l*10*L))> 

K)B  Lj«l  STL^  1  UNTIL  50  no 

WHITE(  LlNf  «FFF»  J*L»FnR     KJ  =  1     e;  T  E  P     1     IJNTI 

lu   un      wl k*i *in. j*l1  » J*L  } ' 

WRMECLlNE.Fa     .J*5l.FU»    Kt»l     STEP     1     UNT 
10    0"       rtlK*I*10. j*5l  ] , j*5l  )  • 
«RI  TE(LINECSMP        l  J  )  » 

Eno! 

ENDJ        ** 

e.no;   *<prIntlIne  ** 

inuaunuuuntduntnuuuumnuunjminunnaju 

REAL    PRUCEUilKE  V(X)j  VALUE       XjRLAL        Xl 

ri  E  (i  I  N 

VIbX*(1*X*(+ALPH0*BI*X>) J 


t     UNTIL    11 


STEP     1 


tuss 

1*4* 
***i 


Eni>1 

BEAL  PRnCE 
HEGIu 


A*************************************************** 

UUHE        Pi)T(X)|   VALUF  Xj   REAL  Xj 


EiMD  ; 


PJT»*X«X*(.5*X*(ALI3*Bia*X)); 


00018200 
0  0  018  3  0  0 
00018400 
00018500 
00018600 
00018700 
00018800 

000l89()0 
00019000 
00019100 
00019200 

no019300 
00019400 
0  0  019  5  0  0 
0001 9600 

00019700 
Q0019H00 
00019900 
00020000 

O002O1O0 
O0020200 
D0020300 
no  020400 

0  0020500 
00020600 
0  0020700 
00020800 
00020^00 
00021 000 
00021 100 
0  0021200 
00021300 
00021400 
00021500 
O0021600 

00021 700 
0  0  0  218  0  0 

00021 900 
0 0022000 
0  0  0  2  210  0 
00022200 

0  0  022300 
00022310 
00022320 
0  0022**00 
0  0022500 
0 0022600 
0  0022700 

***00022800 
0  0022900 

*  1!  t  0  0  0  2  3  0  0  0 
00023100 
00023200 
00023300 
O0023400 

XIIOOO235O0 
00023600 

X  *  X  0  0  0  2  3  7  0  0 
00023800 
0  0  023900 
00024000 
00024100 


IL 
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%%%%%%^y%i>1>t%i>i*i>t^1>ti*%%*%X%t%%**%t%%%%%%%%%*%%%%r%t%%%%%%%t%tt%%%t%%%%0002'i200 

HKOCt'UUNfc-  NYS^(N»Y»Z))  0002*300 

VALUE  N»  0002**00 

INTEGER   M;  0002*500 

REAL  AHKaY  Y,^[0]»  0002*600 

%%%t%%*%1,%*it%*%thl%%%%%%%%%t%%t%%%%%%%%i%%*%%%%%%X%%%%%%X%%%%%%%t%%%%%0002<i700 

t  0002*800 

THIS  PROCEDURE  PERFORMS  ONE  STFP  nF  RUNGE-KUTTA0002*900 


the  equation  is 


METHUqC NYSTRUEMS  FORMULA), 

(0*2)Y/DX*2»niFFUN( Y) 
MITH  INITIAL  v*LUES    Y-ANO  Dy/Dx=Z. 
Y  IS  A  VECTOR  WHOSE  ELEMENTS  ARE 

Y  Co  J . Yll I..... YCnI. 
H  IS  STF.P  SUEtI.Et   Y 

from  y  at  x. 
huth  Ends  are  free. 
un  may  24  1972  k.miura 

K.MIURA 


AT  X*h  IS  EVALUATED 


1     UNTIL    N    00 


*  REVISED 

*  Kt.VlSED    UN    MAY    29 
I 
t, 

H  E  c  I  N 

I P4  T 1. i»  E  R    I  '• 

K5»=K7 »  =Ky  t ,0  * 

n^u«Y[u, 

FUR    I»=l    STEP 
HfCiIN 

Ka»  =  Kl[ni=-^9>(K9;sV(-(YTl«  =  YlPi)  +  {YlPll3Yri 

♦  n ) )); 

KMsKbj  XK6xYUI-l] 

K5l»YH*«4*((Zlll»ZCl3)*.?*K4*H)*H»*l  ST  IT 

K21I  Jl»V(K5-K6) J        X*K?  TEMPORARY  USE  FOR 
*8»=K7>  X  Y2CI-11 

k7I«yIUh2*(ZI1/J*h*k'»/9)  J    X2Nn  ITERATION] 
K3t  I  1  «*V(K7-K8) » 

eno; 
k  b  »  s  o ; 

K2[l  ]  Jsk3[1  j:aK2[N+l ] jsk3[N  +  1] ,sOj 

FUH  I  «  *  1  STEP  1  UNTIL  N  00 
BEGIN 

Y  I  1  *  = 

k2[ T] :=k2[ 1+1 ]-k2[ I ] j 

K3tlJ!»K3[l*l]-KitI]> 

K6«=KS>  %Y300028*00 

K5l«Y[I]*,80*(ZtI]*.20*((?Ml«Kl[Tl)*Y!l)*H)*Hl00028500 

Z[  I  ]  »=**HM23*ZI1/I92j  00028600 

Kl  t  I  J J  =  V(Kb-K6) J       *K1  TEMPORARILY  USED  FOR  V00028700 

LNll  *  00028800 

K 1 C 1 1 !*X1 CN+1 ] :*0;  00028900 

FUR  I:=l  STEp  1  UNTIL  N  00  00029000 

BEGIN  00029100 

KfttsKl  [1*1  ]-Kl  [  I  ]>  00029200 

Y[I)|3*  +  H*(Z[n  +  H*(75*(Yin  =  K2tTn  00029  300 

-(K6tx2^*K3[l))*25*K4)/l9?)|  00029*00 

Zm»=**H*(i25*(YIl+K4)-(K6*K6  +  K6))/i92!         000295oO 

ENDJ  00029600 

YUm-M  J  l»YtN]  j  00029700 

ENOJ    XNYSTROm  FREE  HUUNQARY  REVISED  00029800 

%%%%%%%it*iit%%t%-i,ii.lH%%i%%%i%%t%t%*i>%i%%%%*i%%%%%%%%%%%%%%%%t%i%%%%%%%000299Q0 

%*%%%%%%% t%%*t%**%*t%il%%l%l%%%%%%%%l*%%%%%%%%%%%%%%%%%%%%%%l%%%%%%%%%%XOOO 30000 

00030100 
WRITL(STATI0N.<X10."TYPE  IN  PARAMETERS  IN  THE  FOLLOWING"  00030200 

00030300 


00025000 

00025100 
00025200 

00025300 
00025*00 
00025500 
0002S600 
00025700 
00025800 
00025900 
00026000 
00026100 
00026200 
00026300 
00026*00 
00026500 

00026600 
00026700 
00026800 
00026900 
00027000 
00027100 
VT00027200 

00027300 
00027*00 

00027500 
00027600 
00027700 
00027800 

0002^00 
00028000 
00028100 
00028200 

0002&3QO 
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<  ««     Up     ATJuSVLEfT     TEmpEHATuHL«/«RIGhT     TEmpERAT||RE«/ 

"quadratic  cueff."/"cl)etic  coeff  ,w/"t  i  me  l  i  m  i  t"/"  i  n  i  t  t  al  condition' 
/-     if  randd*  put  i  and  specify  yin  and  Zin  ." 

/"      IF  TkANSMISSIHN  PROPERTIS  ARE  TO  RE  STUDIED  »P(|T  0"/ 
-      IN  i4H'Ch  ^asE  L^-FT  TEmP  ANq  R1&ht  TEmP  ARE  I N I  T  T  A  (.  V^L" 

/"    YIN"/"    ZlN"/"PHlNT  Y&Z(PUT  1  IF  SO  AND  SPEC  I F YsC AL I NGS ) " 

/"    SCALE  Y"/"    SCALE7«/«TIME  I  NCREMENT«/»PR  I  NT  TEMPERATljRF  (  I  F  SO" 

H   SPECIFY  IHE  INTERVAL  a!  wHICH  YUU  "ANT  THE   GRaPH" 

/"BOUNDARY  CHUSS  CHE C  K  I  Nl/V  "H  AZEL  T INE  INTERVALS  If  NO  NEED*" 
/"HAZELTlNF  SCALlN(i(0  IF  NOT  NEE0EU"//>  )  I 

REAn(STATIl)N./.J.AMP.wnTH,ALPHO.RI»TIME,STP.TMP,SlJM.FLG.ISTART. 

ILAST.h.Kl&3.L1»FlG2,hZsCL)> 
WRITE( STATION. <//x10."ThANk  YOU   "/>)! 
wRI TE<LlNE  •<lOR10'3Z^R10*-i>'N»AMP' 

<Jl)TH.TlME»STP»  ALPHU»HI»FLG»  1START.  ILAST» 

n, f-  Li>3.Ll  .Sum,  tmP)  . 

K*XX*****i.CUNSTANTS  i  t,%%t  t  it  Hi  1 1, 1%%1 1  I 

Al  I  3«  =  ALPHu/3» 
K 1 4  «  ■  I?  I  /  4  I 

HM^J  *BM*2  } ; 


HI  «  =1/H» 


HZ »=H«?J 


*  *  *  H     INVERSE 


Twiiib  I  mi**  lb  | 
tiX**Ji***/.  it******4rilTlAL  C  UN  i)I  T  1  DNS     I  l  t t  l  I  it k t%%%%X%%Xt%l% ll%t%% I 
FOR  I««l  STEP  1  lmTIl  1000  ()()  HI3I«RNI1XJ  XHIDLING  (JF 
/MX!t(  AMP*H|)TH)  *5UM» 
IF  STP  E  (j  L  1  THEN  t     THERMAL  SITUATION 

mjr  1 1 »i  step  i  u,nt  II  n  do 

HF.<-,IN 

Y[  I  ]  »  =  TMP*(RNOX-.S) I 


ZMllalF     RND*     GFU     .5     THFN    MAxWFlL(ZMX) 

Else  -mAxwEllc zmx ) j 


End 


ELSE 


*  transmission  property 


HEuIN 

Z[  1  )  '  ,AMp; 

Z[  N  ]  •«-«OTHI 

ENnl 

YC  N*l  1«»Y[ N]  » 

Marti 

FDR  I  I «t  STEP  1  UNT IL  N  DO 
XCU[  I  ]  l»l*6| 

L«  I   IF  n  <  MMt  THEn 

HEbIN 


ooo3oaoo 

00030500 
00030600 

00030700 
00030800 

00030900 
00031000 
00031100 
00031200 

00031300 
00031^00 
00031500 
00031600 
00031700 
00031800 
00031900 
^0032000 
00032100 

0003??00 
O0O3?3OO 
00032400 
00032500 

000  32600 
00032700 

00032800 
00032900 
00033000 

00033100 
00033200 
00033300 
00033400 
00033500 
00033600 

00033700 
00033710 
00033800 
000  3  3900 
0  0  034000 

00034010 
00034100 
00034200 

00034300 
00034400 
00034500 
00034600 
00034700 
00034800 
000349  0  0 
00035000 
00035100 
00035200 
00035300 
00035400 
00035500 
00035600 
00035700 
00035800 
00035900 
00036000 
00036100 
00036200 
00036300 
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COMMENT 
CuMMpNT 


CUmmFnT 


KKJ=K/50+l t 

WRITECLINEISKIP  l3>» 

WHITE(LINE.<X20»"C0LLISIUNS  BETWFEN  TIME»"»I*« 

«   AND  TImEs",  U,//x4."(  LEpT)"  . X4 . " T I mE" , X5  , 

"VELnClTY",X25>'"(RlGHT)'"X6'"TtME','X«i,wvELnCTTY*,/ 
<23,"<FRnM)",XlO»  'Tnl".  X35.*(FRnM>',»Xl0.fT0>',//> 

( K»  K  +  50) j 
rUR  jtxl  STEp  1  uNTk  SO   „0 

h  e  r,  i  n 

Fnrt  JJ»  =  1  STE^  1  UNTIL  HI  DO 
«EGIN 

IF  JJ         MUU   LI  E«l  0  THEN 
REGIN 

IF  Y[i]+H*(ZMXlaZ[l])<..oo01  THEN 

«E(iIN 

LCNTR1»LCNTR*1 ; 

Z[  1  J  J=MAXWp|_L(  AMP)        > 

WRlTE(LlNE,FlF,K*J*H,ZMx,7tl>» 
LFLUXJ=LFLIJX*(ZI1  *ZMX)*(ZT1  -ZMX)» 
ENOJ 
IF  Y[N]+H*(ZMXl*Z[N] >>,0001  THEN 
BEGIN 

RCNTR  JaRCNTH+1 : 

zin  = 

ZIN  J  »="MaX»«ELL(  WDTH)  J 
RFLUXl»RFLUX-(Zn    *ZMX)*(Zll    -ZMX)» 


C  U  M  M  FN  T 


WRlTE(LlNE»FRF,K  +  J*H#ZMX,Zini 


ENnj 


ENn    Uf    COLLISIONS     *ITH    H^AT    RESErv^^I 
■JYSalN        »Y»Z>J 
FOR     I:=l     STEP     1     UNTIL    N    DO 
KNTp[ I ] :=*♦/[ I]**2; 

end  of  jj  loopj 
if  flg  eol  1  then   'print  y  and  z 

begin 

Y(  N-f  1  ]  tsY[N]  ;  SOMjsO; 
ZMX J=0{ 

FOR  1»  =  1  STEP  1  UNTIL  N    DO 
BEGIN 


*ii  «»tQ[  Jiii  *  =  zr  1 1 » 
Sum;s**,5*(  zn*zi  UZmx 
+  czmx »=poT(TPr j»  n i 

Yl  1  +  1  J-YC  I  1  ))  )  J 
ENOJ 
TP[J.NJI»Y[1J> 
TPC  J»N+1 3  JcYtN] X 

TQ[ J.N+1 1 JsSUMJ   X  TOTAL  ENFRGY 
END  Of  SAVING  DATA  rOR  PRlNTlNr.J 
IF  FLG2  NE«  0  ™EN  I F  J  M0D  FLr>2  F0L 
BEGIN   XX  HAZFLTINE 

FOR  I»=l  STEP  1  UNTIL  N  On 

YCOt  m(Z[I  ]*HZSCL> 
HAZELPLUT(XCO.YCO»N,K*J); 
END  HAZELTINEJ 
END  OF  J  LOOPJ 


0  THEN 


00036400 
00036500 

00036600 
00036700 
00036800 
00036900 

0003r000 
00037100 
00037200 
00037300 

00037400 
00037500 
00037600 
00037700 
00037800 
00037900 
00038000 
0003*100 
00038200 
00038300 
00038400 

00038500 
00038600 

00038700 
00038800 
O0038900 

00039000 
00039100 
00039200 
00039300 
00039400 
00039500 
00039600 
00039700 

00039800 
00039900 
00040000 

000*0100 
00040200 
00040300 
00040400 
00040500 
00040600 
00040700 
00040000 
00040900 
00041000 
00041 100 
00041200 
00041300 
00041400 
00041500 
00041600 
00041700 
00041800 
000«1900 
00042000 
00042100 
00042200 
00042300 

000*2^00 
00042500 
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LF[KK j »=LFLUX; 
«FIKK  ]  JsRfLUX J 
If  FLG3  NEQ  0  THEN  IF  K  MOD  FLG3  EQL  0  THEM 
KINFTICTEHP(N.KNTP.K*50  .N){ 
IF  FLG  EQL  1  THEN    SPNINT  Y  AND  Z 
»  E  G  T  N 

polygraphs, n  .tp.istart)!   %     coordinate 
polygraphs, n  . t« . i l ast ) >  **velnctty 

eno» 

IF  (jl=K+5(>)  MOO  1000  EQL  0  THEN 

hEg  iN 

WHI TLCLINF (SKIP  i 1 ) I 

w*ITt(LlNE.<MO,"FLUx  LEAST  SQUARE  FITTING  AT  ".15/ 
X5."LEFT«.x25."RIGHT'V>.J)J 

pl)K  L,=100  STEp  100  UNTIL     J-500   DO 

HEUl  N 

LSStHLF.L.  J  ,50»SUM,ZMX)  J 

LSS^(RF.L.J  ,b0.SUM,Zi i ) • 


"KITE  (I.  INE.<2(Xb»I5.X2.Ell.a)>»L.ZMX*tS, 
L*2ll*t 5}  J 

END  (tF  FLUX  LEAST  LINLaW  F  I  T  T  I  N  G  J 

w»lTE(STATIUN.<//,,MuHl'?    I  F  SO.  TYPE  1"///>)I 
KEAD(STA1 IUN./.M) J 
lp  M  NEy  1  THEM 
GO  TU  Lb» 

tNDJ  ****  MOO  1000 
K  tiK  +  Su1 
GU  TO  L4; 
EN()|gft*s*  lip  k  L00P(50  SEC)) 
L5» 

END. 


00042600 
00042700 

0004?800 
00042900 
00043000 
00043100 
00043200 
00043300 
00043400 
00043500 

00043600 
00043700 

00043^00 
00043900 
00044000 
000^4100 
00044200 

00044300 
00044400 
00044600 
00044700 
00044*00 
0  0044900 
00045000 
00045010 
00  045  100 
00045200 
00045300 
00045400 
00045500 
00045600 
00045700 

00045600 
00045900 
00046000 
00046100 
00046200 
00046300 
00046400 
00046500 
00046600 
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